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FOREWORD 

MAPSEP (Mission Analysis Program for Solar Electric Propulsion) 
is a computer program developed by Martin Marietta Aerospace, Denver 
Division, for the NASA Marshall Space Flight Center under Contract 
NAS8-29666. MAPSEP contains the basic modes: TOPSEP (trajectory 

generation), GODSEP (linear error analysis) and SIMSEP (simulation). 
These modes and their various options give the user sufficient flexi- 
bility to analyze any low thrust mission with respect to trajectory 
performance, guidance and navigation, and to provide meaningful sys- 
tem related requirements for the purpose of vehicle design. 

This volume is the first of three and contains the analytical or 
functional description of MAPSEP. Subsequent volumes relate to opera- 
tional usage and to program logical flow. 


! 

i ' 


4 



f 



v-> 


Acknowledgement 


The authors wish to express their thanks to Jeanette Wearden for her 
skill and stamina in typing all three volumes, to Susie Borah for her 
art work on the flow charts, and to Kathy, Paula, Joan and Marion for 
their patience and understanding. 


TABLE OF CONTENTS 


Foreword 
Acknowledgement 
Table of Contents 

1.0 INTRODUCTION 

2.0 PROGRAM DESCRIPTION 

3.0 NOMENCLATURE 

4.0 TRAJ 

4.1 Equations of Motion 

4.2 Trajectory Termination 

4.3 Trajectory Accuracy 

4.4 Trajectory Repeatability 

4.5 State Transition Matrix Generation 

4.6 Covariance Integration 

5.0 TRAJECTORY GENERATION - T0PSEP 

5.1 Nominal Trajectory Propagation 

5.2 Grid Generation 

5.3 Targeting and Optimization 

6.0 LINEAR ERROR ANALYSIS - G0DSEP 

6.1 Augmented State 

6.2 Covariance Propagation 

6.3 Measurement Types 

6.4 Filter 

6.5 Generalized Covariance 


Page 

ii 

iii 


iv 

1 


4 

10 

13 

13 

22 

23 

25 

27 


30 


34 


35 

35 

36 

55 

56 
58 
63 
76 
80 


6.6 


Guidance 


82 





7.0 TRAJECTORY SIMULATE 'N - SIMSEP 90 

7.1 Program Scope and Methods 90 

7.2 Definitions and Concepts 93 

7.3 Guidance 95 

7.3.1 Linear Guidance 96 

7.3.2 Linear Impulsive Guidance 97 

7.3.3 Low Thrust Linear Guidance 100 

7.3.4 Non-Linear Guidance 103 

7.4 Simulated Orbit Determination 108 

7.5 Thrust Process Noise 109 

7.6 Guidance Execution Errors 110 

8.0 REFERENCES 112 

9.0 APPENDICES 113 

9.1 Conic Equations for Position and Velocity in 

Elliptical and Hyperbolic Orbits (Appendix 1) 113 

9.2 A Generalized 4th Order Runge-Kutta Algorithm 

with Range's Coefficients (Appendix 2) 115 

9.3 Newton's 3rd Order Divided Difference Inter- 
polation Polynomial (Appendix 3) l 1 7 

9.4 Analytical Expressions for Terms in F Matrix (App. 4) 1.9 

9.5 TOPSEP Injection Modeling (Appendix 5) 124 

9.5.1 Injection Controls 124 

9.5.2 Tug Multiple -Impulse Orbit Transfer 129 

9.6 Control Weighting Matrices (Arpendix 6) 137 


9.7 Integrated State Transition Matrices for Computing 
the Targeting Sensitivity Matrix (Appendix 7) 




1. INTRODPCTION 


A major requirement for spacecraft systems design is the effec- 
tive analysis of performance errors and their impact on mission success. 
This requirement is especially necessary for low thrust missions where 
thrust errors dominate all spacecraft error sources. Fast, accurate para- 
metric error analyses can only be performed by a computer program which 
is efficiently constructed, easy to use, flexible, and contains model- 
ing of all pertinent spacecraft and environmental processes. MAPSEP 
(Mission Analysis Program for Solar Electric Propulsion) is designed 
to meet these characteristics. It is intended to provide rapid evalu- 
ation of guidance, navigation and performance requirements to the 
degree necessary for spacecraft and mission design. 

The baseline design of- MAPSEP was taken from a previous study 
effort (Reference 1). Suitable modifications to the design were made 
which reflected subsequent operational experience (References 2 and 3) 
and actual construction and testing of MAPSEP. Considerable knowledge 
was also gained in the construction and usage of the engineering ver- 
sion of MAPSEP which was actually three separate programs that corre- 
sponded to the modes (Figure 1-1): TOPSEP, GODSEP and SIMSEP. Driving 
considerations in the program design and construction were: realistic 

vehicle and environment modeling consistent with preliminary vehicle 
design, flexibility in usage, computational speed and accuracy, mini- 
mum core utilization (for turnaround time and operating costs) and 
maximum growth potential through modularity. 




Figure 1*1 MAPSEP Modes 
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This document is the first of three volumes. Contained herein 
is a brief introduction to MAPSEP organization and detailed analyti- 
cal descriptions of all models and algorithms. These include, for 
example, trajectory and error covariance propagation methods, orbit 
determination processes, thrust modeling and trajectory correction 
(guidance) schemes. This analytic background is necessary to fully 
understand program operation and to maximize program capability with 
respect to the user. 

The second volume is a description of program usage, that is, 
input, output, recommended operating procedures and sample cases. 

The third volume is a detailed description of internal MAPSEP struc- 
ture including macrologic, variable definition, subroutines and logi- 
cal flow. 


2. PROGRAM DESCRIPTION 


This section summarizes MAPSEP' s function and use, and MAPSEP 
structure. These areas are discussed in greater detail in the user's 
manual and programmer's manual, respectively. 

As mentioned earlier, MAPSEP is composed of three primary modes. 
Each mode is intended to serve a given function in the mission design 
sequence. TOPSEP (Targeting and Optimization for SEP) is used to 
generate numerically integrated trajectories consistent with dynamic 
and system constraints. Performance data and related sensitivities 
are computed in the process of trajectory generation but can also be 
obtained by parametric application of TOPSEP. Indeed, each mode 
readily lends itself to parametric use which is a necessary feature 
for mission and system design studies. 

GODSEP (Guidance and Orbit Determination for SEP) is used to 
perform a linear covariance analysis about a selected reference tra- 
jectory, generated by TOPSEP. Various dynamic and measurement related 
error sources are applied in a statistical sense to compute trajectory 
error covariances. These covariances, corresponding to estimation un- 
certainty (knowledge) and to actual trajectory deviations from the 
nominal (control), are propagated through a sequence of mission events 
thruster switching, navigation measurement /state update, trajectory 
correction (guidance), etc. Thus, GODSEP computes a time history of 
the ensemble of all expected trajectory errors, and in the process 
displays such useful system parameters as required thrust control 
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authority, predicted terminal miss, additional fuel expenditures for 
off-nominal performance, etc. 

SIMSEP (Trajectory SIMulation of SEP) is used in the latter 
stages of system design. It deterministically simulates the trajec- 
tory including the application of discrete dynamic errors. Trajectory 
corrections are simulated in an operational sense through a thrust 
update design and execution process. Navigation is simulated by 
sampling estimation error covariances (generated by GODSEP) prior to 
each guidance event. By operating SIMSEP in a Monte Carlo fashion, 
any desired number of simulated missions can be obtained, and estimated 
singly or collectively in statistical displays. 

A fourth "mode", REFSEP (REEerence SEP), is actually an t3xpansion 
of TOPSEP to provide a greater amount of trajectory and navigation 
related data for a particular reference trajectory. 

Each mode of MAPSEP uses a commc trajectory propagation routine, 

TRAJ. This guarantees trajectory reproducibility among modes. TRAJ 
integrates the equations of motion in Encke form using a fourth order 
Runge-Kutta scheme. Covar 4 ance propagation and transition matrices 
are computed by integrating variational equations simultaneously with 
the equations of motion. An option exists in GODSEP which stores a 
complete set of trajectory parameters and transition matrices as it 
is generated by TRAJ onto a magnetic disc called the STM file so that 
subsequent error analyses will not have tc regenerate the data. The 
information on disc can then be transferred to magnetic tape for per- 
manent storage, if desired. A more limited option is available for 
TOPSEP and SIMS^ which stores only the initial (input) trajectory 

data on disc. ; 

/ 
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Input to MAPSEP Is primarily through cards using the NAME2 CY 
feature, with supplementary means depending upon mode and function 
(Table 2-1). All modes require the $TRAJ namelist which defines 


Mode 

INPUT 

OUTPUT 

Namelist Formated Tape 

Cards (or disc) 

Punched Tape 

Cards (or disc) 

TOPSEP 

$TRAJ None STM 

$T0PSEP 

None STM 

GAIN 

i— . . . .... 

GODSEP 

$TRAJ 

$G0DSEP Event STM 

SGEVENT Data GAIN 

States 

Covariances STM 

Guidance GAIN 

SUMARY 

SIMSEP 

$TRAJ None STM 

ISIMSEP 

0GUID 

Statistics STM 

GAIN 
SUMARY 

REFSEP 

$TRAJ Print STM 

Events 

None STM 


TABLE 2-1. MAPSEP Input /Output 

the nominal trajectory and subsequent mode usage. However, If 
recycling or case stacking is performed It is not necessary to in- 
put $TRAJ again unless desired. The second namelist required for 
each mode corresponds to mode peculiar input and bears the name of 
that particular mode. Additional namelist, formated cards, and tape 
input are generally optional. Besides the standard printout asso- 
ciated with MAPSEP auxiliary output can be obtained which will 


facilitate subsequent runs 
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The structure of MAPSEP is organized into three levels of 
"overlays" which are designed to minimize toi :1 computer storage. 

At any given time, only those routines which are in active use are 
loaded into the working core of the computer* The main overlay 
(Figure 2-1) is always in core and contains the main executive, 
MAPSEP, and all utility routines that are common to the three modes. 
The primary overlays contain key operating routines of each mode, 
that is, those routines which are always needed when that particular 
mode is in use. Also included as a primary overlay Is the data 
initialization routine, DATAM, where $TRAJ namelist is read, trajec- 
tory and preliminary mode parameters are initialized, and appropriate 
parameters are printed out. 

The secondary overlays contain routines which perform vari< 
computations during a particular operational sequence. Included -re 
data initialization routines, analgous to DATAM, which operate on 
mode peculiar input and perform mode initialization. An example of 
core usage in the changing overlay structure may be provided by a 
standard error. analysis event sequence. Error analysis initializa- 
tion is performed by the overlay DATAC. Transition matrices are then 
read from the STM file, the state covariance is propagated to a 
measurement event, and the overlaj MEAS is called, which physically 
replaces, or overlays, the same cort used previously by DATAC.. 
Similarly at a guidance event, overlay TRAJ will replace MEAS to 
compute target sensitivity matrices and ovtrl^y GUID will then 


> 







Figure 2-1. OVERLAY STRUCTURE 
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replace T RAJ to compute guidance corrections. Overlay switching is 
performed internally and is transparent to the user. 

All of the routines and structure of MAPSEP are constructed to 
minimize core storage (thus reducing turn-around time and computer 
run cost) yet retain the flexibility needed for broad analysis re- 
quirements. Furthermore, routines are built as mouular as possible 
to reduce the difficulties in future modifications and extensions. 
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3. NOMENCLATURE 

The following symbols are used throughout the Analytic Manual, and 
to a great extent in the User's and Program Manuals. However, deviations 
from these symbols may occur in localized discussion if required for pur* 
poses of clarity. 


SYMBOL 


DEFINITION 


propulsive acceleration 

propulsive exhaust velocity 

cross covariance between i and j parameters 

target error index 

dynamic variation matrix 

performance gradient or thrust transfor- 
mation matrix 

observation sensitivity matrix (WRT state) 
Identity matrix 
filter gain matrix 
spacecraft mass 

covariance of state deviations or elec- 
trical power or projection operator 

effective power at 1 AU 

dynamic (thrust) noise matrix 

spacecraft position 

solve-for parameters 


target sensitivity matrix WRT control para- 
meters 





n 


DEFINITION. 

event time, or target variables, or 
thrust 

dynamic (consider) parameters 

thrust acceleration proportionality (throttlin 
Control parameters 

spacecraft velocity 0£ measurement para- 
meters 

ignore parameters 

weighting matrix 

spacecraft state 

guidance matrix 

propulsive efficiency 

transition matrix of dynamic parameters 

gravitational constant 

relative range 

standard deviation 

correlation time of thrust noise 

transition matrix of augmented state 

state transition matrix from time t. to 

t k 

k+1 

thrust noise 

DEFINITION 

state control covariance 


matrix evaluated over time interval t 
to Stfi 


k 


state knowledge covariance 



SUBSCRIPT 

DEFINITION 

( >» 

planet related parameters 

< >, 

solve- for parameters 

< >» 

measurement consider parameters 

< >. 

ignore parameters 

( > x 

spacecraft state parameters 

MISCELLANEOUS 

DEFINITION 

G&N 

guidance and navigation 

OD 

orbit determination 

PGM 

Projected Gradient Method 

S/C 

spacecraft 

SEP 

solar electric propulsion 

WRT 

' with respect to 

E[ ] 

expected value operation 

( ) + 

post-event value 

( )* 

pre-event value 

• 

< > 

time derivative 

< A > 

unit vector 
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4.0 TRAJ 

Essential to any program used for performance and navigation 
analysis is an accurate, but computationally efficient, trajectory 
propagation routine. This routine must contain realistic models of 
the dynamic processes acting on and performed by the spacecraft. In 
MAPSEP, the subroutine TRAJ fulfills this role. TRAJ is designed to be 
used by the three modes TOPSEP, GODSEP and SIMSEP, and is capable 
of duplicating the same trajectory in all modes. 

The trajectory overlay TRAJ propagates planetary and inter- 
planetary low thrust trajectories, using Encke's formulation of 
the equations of motion, from any epoch to a termination condition. 

TRAJ can optionally propagate the state covariance or the state 
transition matrix along the trajectory for either the basic state 
r. an augmented state. Two of the most important features incorpo- 
rated into TRAJ are the variable integration step algorithm and 
trajectory repea tibility. 

4.1 Equations of Motion 

In Encke's formulation, all accelerations other than those due 
to the gravity of a primary body are called perturbing accelerations. 
Position (r^) and velocity (r^) vectors are computed relative to the 
primary body using two body formula. The deviation vectors from 
the reference conic position and velocity vectors, $ r_ and 5 £» 
respectively, are the direct results of numerically integrating the 
sum of the perturbing accelerations. The true position and velocity 
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vectors are, respectively, 

r * r + Sr 

“ -c — 

® • . 

r = r + % r 

— — c — 

Since r^ and r^ can easily be computed from conic formula (See 

* 

Appendix O, the only problem is to compute £ and £. 

• * 

Let &£ be the acceleration deviation from two body motion, such 

• « 

that &r is the sum of all perturbing accelerations, e.g., other 
bodies, thrust, etc. 

N 

t+ ^ [r ♦f(« 1 ) % ] 


+ i + % 


The first term is the difference between two body and perturbed two 
body motion. The second term is the sum of the accelerations due to 
the N perturbing bodies. The third term (£) is the acceleration due 
to thrust. The fourth term (a ) is the acceleration due to radiation 
pressure. 

The first term is computed from the following equations (See 
Reference 4), 

1 + <1 + ot) ' 


U , (St " 2r) • £ £ 

2 

r 
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where r is the true s/c position vector relative to the primary body 
and ja is the gravitational constant of the primary body. 

To compute the second term, the heliocentric position vectors, 
r^, of the perturbing bodies, are computed from mean analytical 
orbital elements to obtain 


gi * r 


r + r 


r . 




°^i « 


f, 


L 1 + (1 + 0^) J/z J 

f L . 2 I • Pi 

L '* 


where £p is the heliocentric position vector of the primary body. 
The acceleration vector due to radiation pressure is 


%L 


1.024 x 10 8 A C r 
2 

© r 


A 

r 


where 1.024 X 10 8 


r 

A 


The option exists 


- solar flux constant 

• instantaneous mass of the s/c 

- heliocentric position vector of the s/c 

• effective cross sectional area of the s/c 

- coeffient of reflectivity. 

in TRAJ, to include or exclude the effects of 
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radiation pressure when propagating both planetary and interplanetary 
trajectories. Before defining the acceleration due to. thrust j|, we 
will model the power subsystem. There are two power subsystem models 
used by TRAJ for low thrust. They are solar electric and nuclear electric. 
The power to the thrusters (p) is 



£- p L (t-t DL )J 



solar 

electric 


P 




P 

max. 


if P 


> 


P 

max 


or 


< r 


min 


solar 

electric 


exP [* P i J ( t“ t DL ) ] - P HK 


nuclear 

electric 


P - Power available (at 1 AU for solar, at energization for 
° nuclear electric 

- (Empirical) Constants defining solar array characteristics 

r - Heliocentric position magnitude of the S/C 
P^ - Power decay constant 

t - Time from epoch 

- Time delay 

Pgg - Housekeeping power 

P - Maximum allowable solar electric power 

max 

r , - Heliocentric distance at which P reaches P 

sin max 

The exponential term in the solar electric expression describes the 

degradation of the solar array as a function of time. 


! 

i 
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The thrusters provide the spacecraft with the ability to maneu- 
ver. By changing the orientation of the thrust vector, and maintain- 
ing that orientation for a long enough time, it is possible to steer 
the spacecraft and "shape" the trajectory. The thrust controls are 
defined in terms of constant parameters over a given time segment of 
the trajectory. Thus the user or mode can specify the following 
controls for each segment (up to 20) : 

1. Thrust policy: Cone + Clock (sr. .»rientation) , In and 

Out of Plane Angles (orbit plane coordinates), or coast. 

2. Segment end time (referenced to launch) of the current 
segment. 

3. Throttling level, 

4. Cone angle (or In Plane Angle). 

5* Clock Angle (or Out of Plane Angle). 

6. Time rate of (4). 

7. Time rate of (5). 

8. Number of operating thrusters. 

The thrust policy is merely thrusting or not thrusting (coasting). 
During thrust, the user has an option as to the reference system for 
the acceleration vector. The two reference systems are Cone and Clock 
Angles, Figure (4-1), and In and Out of Plane Angles, Figure (4-2). 

The latter is primarily used for Earth orbital applications and the 
former is primarily used for Interplanetary missions. 


4 
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(-SIM) 



The spacecraft is assumed 
to be oriented with the Z' 
axis in the sun-spacecraft 
line and the X' axis in a 
plane containing Z' and Z 
(reference star directionf. 


Figure (6-1). Cone and Clock Angles 



6 - Out of plane angle 
T - In plane angle 


Figure (4-2). In and Out of Plane Angles 


/ 
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The acceleration due to thrust is 

|.|- °-° 02 t F “T 

tn c 

.002 - Conversion factor 

- Averaged efficiency of the thrusters 

P - Power delivered to the thrusters 

U - Throttling level 
T 

n - Spacecraft mass 
c - Exhaust velocity 

The mass is numerically integrated from the equation 

* - ma 

m * 

c 

The thrust acceleration vector a/relative to the spacecraft for 
the two reference systems is 

• ' « a cos CLOCK sin CONE 

x 

a' y - a sin CLOCK sin CONE 

a' « a cos CONE 

m 

for the Cone-Clock system and 

a* x • a cos S cos V 

• 'y * a cos $ sin T 


a sin S 



I 


t 



1 
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for the In and Out of Plane system. In order to transform a/ into a, 
which is in heliocentric ecliptic coordinates, we must, define the 
matrix A such that 


where 


a ■ A a* 


i 

O 1 


o 


a - £x’ : y ‘ ■ r ' j 

and V - r/|r| 

Y' ■ rxZ/irxZ i 

— g i — —5 i 

X' ■ x Z ' 

jr is the heliocentric spacecraft position vector and Z are the 

“"S 

ecliptic direction cosines of a reference star for the Cone/Clock 
System. For the In/Out of Plane system, A is defined in terms 
of the position and velocity vectors relative to the primary body, 
let £ be the position vector and v be the velocity vector, then 
A can be defined as 


where 


and 


* -[*! i i 2 j A 3 1 


“ V / |v 1 


-3 


r x v 
lr x v | 


A- A 
“2 -3 


x A j 


„ - j 


A ^ 
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Since the trajectory is made up of segments, TRAJ propagates 
the trajectory using the appropriate set of controls. over the inter- 
val the controls are in effect. Updates are automatically performed, 
•t the beginning of each new segment. 

• • 

Mow that the perturbing acceleration, fcr, , is defined, we can 

• • • 

obtain Sr. and Sr . To do this, £ r, is numerical' tegrated 
with a generalized 4th Order Runge-Kutta algorithm f • . rst order 
differential equations (Appendix 2). We can express as a set 

of first order equations 


• Si 

Si ■ 


%x can be numerically integrated to give £x(t) 
given the following initial conditions: 


■m 


when 


at t ■ t 


>• Ss ■[!] 


* nd ^ ■ to . to ■ t. 


The propagation of $£. (and Sr ) continues until Sr, is greater 

than or equal to some prescribed value Sr , then "rectification" 

max 


occurs. 


Hence, when £r ^$r . at some time, t, 

max 

• • 

ve reinitialize r • r(t), r * r(t) and S 




and compute a new reference conic orbit. Rectification ensures that 


i b 


the conic is always "c'ose" to the true state. The propagation 


i 
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continues until Sr is again greater than or equal to Sr . and 

max 

so on to the end of the trajectory. 

The use of constant thrust controls over discrete time intervals 
stakes the trajectory discontinuous in acceleration at the control 
switching boundaries. The Integration is thus done piecewise. That 
is, the state at the boundary between segments is used as initial 
conoitionc (rectification) for the next segment. 

4.2 Trajectory Termination 

There are four options for terminating trajectories in TRAJ: 

(1) final time, (2) closest approach to a target body, t3) sphere of 
influence of the target body, and (4) a radius relative to a target 
body. Termination at final time is straight forward. For the other 
conditions, termination criteria are tested after each integration 
step. Once a cutoff condition is sensed, a step-size is computed so 
that TRAJ can propagate *.o an interpolated time. TRAJ takes the three 
previous planet relative position magnitudes plus the present relative 
position magnitude, and the corresponding trajectory times, and fits 
S third order polynominal through the four data points, using Newton's 
3rd order divided difference interpolation polynomial (Appendix 3). 

The independent variables for sphere of influence and stopping radius 
termination are the position magnitudes and the corresponding trajec- 
tory times are the dependent' variables. Since the radius of the sphere 
of Influence or the stopping radius is known before hand, the infor- 
Mt ion that is needed is the time at these position magnitudes. For 






J 

— - I 


c 


o 


u 
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closest approach, the same information is stored as before, but now 
the independent variables are the trajectory times and the position 
magnitudes are the dependent variables. Instead of knowing a time 
for which the position magnitude is a minimum, a value is computed 
corresponding to the minimum of a 3rd order polynomial. 

Of particular note is stopping on the sphere of influence or 
radius of closest approach (Figure 4-3). These stopping conditions 
are used to evaluate B- plane (impact plane) parameters, B»T and 
B*R (Figure 4-4). These parameters form a convenient set of variables 
for the description of the approach geometry for interplanetary mis- 
sions. Let V_ denote the hyperbolic excess velocity of the space- 
—he 

craft. Then 


S 



S X k 

- “ | S X k | 

R = S X T 

where k is a unit vector normal to the reference plane, usually the 
planetocentric ecliptic plane. 

4.3 Trajectory Accuracy 

As with all problems that require numerical integration, some 
criteria must be used in determining the nominal integration stepsize. 
The stepsize algorithm used in TRAJ is empirical, and it meets the 


i 
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requirements of reasonable numerical results and computer run time. 
The algorithm is •: 

h « 6 • |f| % 

where h is the integration stepsize, f is the gravity gradient (See 
Appendix 4) and €. is a user input scale factor. There are con- 
straints on h such that 

h £ 5 days 

for heliocentric trajectories and 

h < 1 day 

for planetary trajectories. 

The effects of h on the state transition matrix are small com- 
pared to the spacecraft s^ate. Mass and mass variation are also not 
strongly affected by h. because they are affected primarily by the 
spacecraft heliocentric position. Therefore, a good choice of h 
ensures a satisfactory trajectory. 

4,4 Trajectory Repeatability 

A major goal in building TRAJ was the ability to reproduce the 
same trajectory in all modes, given the same initial and spacecraft 
related data. To accomplish this, all propagation that resulted in 
altering the nominal integration stepsize was decoupled from the 
event stepsize logic. That is, event times and print times do not 
alter the norminal stepsize; instead, the information at the previous 
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nominal stepsize, t, is saved and a separate stepsize is computed 
for the event or print time, t + (Figure 4-5). If there are 
many events or prints between t and t + h a stepsize h £ is computed 
for each. Afterward, TRAJ returns to the nominal trajectory and 
continues the propagation. This suggests that, if the trajectory 
starts at any nominal integration step the trajectory will be 
duplicated, especially with respect to terminal conditions, for any 
run with different print and event times. For trajectories that do 
not start at a nominal integration step, there will be slight devi- 
ations in the terminal conditions from the nominal, but they will be 
close. 



A 

h 


nominal integration time 

event or print time 

evert or print integration step 

nominal integration step 


Figure 4-5 


Trajectory Preservation 
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4 . 5 State Transition Matrix Generation 

In a linear analysis, the state transition matrix, y» 

is used to map perturbations about tne reference trajectory from one 
epoch to another according to the equation, 


k+1 


fk+1’ k 


* 2 * 


Because this mathematical operation is repeated many times, it is 
important to use the most efficient and accurate method of computing 
<§ available. This is best done by simultaneously integrating the 
variational differential equations which generate J and the S/C 
equations of motion. In MAPSEP these variational equations have been 
implemented and { is augmented to the state variabl s in the inte- 
grator. 

The origin and composition of the variational equations are best 
understood by first considering the equations of S/C motion written as 
a system of first order, coupled differential equations; namely, 

x = f (x, t) (4-1) 


In 4-1, x is a six-dimensional state vector of position and velocity 
components, and f is a vector function giving the time derivative of 
each state component. The most general form for £ may be written as 


f 


v 

& + £L + + f 


(6x1) 


where & is the gravitational acceleration due to the primary body, 



» 

f 
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gf is the gravitational acceleration contribution from secondary bodies, 

and a^ is the thrust caused acceleration term discussed in Section 4.1. 
-T 

^ is a term of miscellaneous accelerations due to radiation pressure, 
planetary oblateness, etc., and is normally neglected for interplanetary 
missions as a low order effect. It should be noted that the three 
primary terms, £, £ and a^, are all dependent on the S/C position vector. 
For example, £ and & depend on position through the law of gravity; a^, 
through the electric power function and the transformation which relates 
the body axis syt.cem to the inertial representation. 

To derive the variational equations, the vector function, £, is 
expanded in a Taylor series about some reference solution to equation 
4-1. Hence, the right hand side of 4-1 becomes 

L (x + £ x> t) = f (x, t) + ^ ~ $ X + & ( $x 2 ) 

i x 

where 6x is relative to the reference trajectory state at time t. 

By neglecting second and higher order terms, this expression reduces to 

Sx - F Sx (4-2) 

where S x is defined by 

Sx = f (x + S x» O - 1 (& t) 


and F is a matrix of first order partials of the vector function £ with 
respect to state components. Explicitly writing the F matrix, it is 
seen that it has only two non-zero partitions. 


F 


33 , 

1 

1 

. hi 

33 

1 

I 

°33 


(6x6) 
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where is a 3x3 identity matrix and is a 3x3 matrix of time 
varying expressions evaluated along the reference trajectory. The 
components of f^, in terms of state variables, are obtained by 
analytically differentiating g* and a^ with respect to r. 

The solution of 4-2 is known from the theory of linear dif- 
ferential equation® to be of the form 

S X = I (4-3) 

where $ is identified as the state transition matrix and is 
representable as 


$ - »■ y. z.v x . y y 

V V V V,' V, v .o> (6x6) 

As noted before, § jc^ is a perturbation, or deviation, from the 
reference trajectory at the initial epoch and, as such, it is arbitrary 
but constant. Differenting 4-3 with respect to time and substituting 
the resulting expression into 4-2 yields 
♦ 

$ - F } (4-4) 

This equation is the variational differential equation for the state 
transition matrix and is numerically integrated to obtain $ over an 
arbitrary interval, t Q to t^, with initial conditions at t Q being given 
by 


\ - 


I 


66 ’ 


From the more general point of view, equation 4-1 can be 
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considered to be dependent on other parameters as well as the usual 
state variables. For example, the trajectory generated as a solution 
to 4-1 is dependent on the thrust controls, the inertial states and 
gravitational constants of planetary bodies, and the solar gravitational 
constant* When some of these parameters are of interest in a linear 
analysis, the state vector .of dynamic parameters is said to be augmented, 
thus increasing its dimension to as great as seventeen. This occurs 
when there are three thrust controls, six ephemeris elements, and two 
gravitational constants in addition to the vehicle state. 

For the sake of modeling simplicity, MAPSEP allows ephemeris 
elements for only one planet to be augmented as dynamic parameters. 

This body is referred to as the "ephemeris body" and is usually 
selected to be the planet that most strongly influences the S/C tra- 
jectory (other than the earth). In many cases, the ephemeris body is 
the same as the target planet. 

For the problem where the state vector is augmented with para- 
meters, the equations of motion are more suitably written as 

_L a - I A «, Xp> /X p. JK 8 ; t) (4-5) 

where x ^ is the augmented state vector, i.e., 
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In the above expression, x corresponds, as befote, to the S/C state 
vector; u , to the thrust control vector; x p , to the ephesteris 
planet state vector; and the jx ' s refer to gravitational constants. 
Following the previous analysis, it is possible to expand 4-5 to obtain 
variational equations for the augmented state transition matrix. This 


differential equation is 

of the 

form 




*A ■ 

f a 

»* 




(4-6) 

where is partitioned 

as 
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and where F^ is given as 
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Explicit definitions for the terms in are given in Appendix 4. Terms 
appearing in are explicitly defined as follows: 


(partials of state component w.r.t. state components) 


(partials of state components w.r.t, thrust controls) 


% * (partials of state components w.r.t. ephemeris planet state 

^ components) ■ 



M * (partials of state components w.r.t. ephemeris planet gravita tional 
^ constant) * 



K * (partials of state components w.r.t. solar gravitational constant) * 


a. 

^Ms 61 


f 



(partials of ephemeris planet state components at the final epoch 
w.r.t. ephemeris planet state components at the initial epoch) » 
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66 
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M' * (partials of ephemerls planet state components w.r.t. ephemerls 
planet gravitational constant) = 


61 


and M' * (partials of ephetneris planet state components w.r.t. so_jr 
gravitational constant) * 

£ x 

61 

Before concluding this section, it should be noted that MAPSEP 
has program logic which allows arbitrary augmentation of dynamic para- 
meters. That is, the program organizes and integrates matrices in 
equation 4-6 dimensioned to accommodate only those parameters requested 
during input. In this way, there is no time wasted in unnecessary 


calculations . 



4,6 Covariance Integration 


In any linear error analysis, a major problem is the prop- 
agation of state error covariances (P) from one event to the 
next event. Two methods are generally used: propagation with 

state transition matrices and numerical integration of the 
covariance matrix differential equations; both of which are 
options in TRAJ. In the latter case, the covariance is inte- 
grated to an event, where operations are performed on P, and 
the updated P is integrated to the next event. 

Given the nonlinear equations of motion 

x - x (x, u, t*) (4-7) 


where x is the spacecraft position and velocity, u are constant 
spacecraft controls and 6j are time-varying thrust parameters 
(nominally zero), these equations can be linearized about a 
reference trajectory such that 
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where £x, §x, 5y. and are errors In the respective dynamic para- 
meters. Both $u and £,u.'are described in terms of the 3x1 parameter 
set: u , cone and clock. The 6x1 £(£ actually models two distinct 

processes for each parameter set (See also Page 62-A). Whereas Equation 
4-7 describes motion of the deterministic reference trajectory, Equation 
4-8 describes the linearized propagation of trajectory deviations result- 
ing from dynamic and a. priori uncertainties. The covariance Integrated 
by TRAJ not only maps dynamic errors but also measurement related 
errors, specifically, uncertainties in three station locations. The 
augmented state covariance is defined as 

p = e [S*a Sil ] 


where 



Sr 

$v 

Sh 

Sii 



so that P - FP + PF T + Q 

where F is similar to that used in the definition of the state transi- 
tion matrix, and is evaluated along the reference trajectory, and Q is a 
process noise matrix. The augmented F matrix is defined as 



1 


I -1 


F ■ 


where I is a 3 x 3 identity matrix, 


ui*i 


and h is the matrix of process noise correlation times 


Analytical equations for terms in the F matrix appear in Appendix 4, 


The process noise, Q, is modeled as a stationary first order Gauss- 


Harkov process and is defined as 



t 

i 
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For a discussion of Q, see Chapter 6 (GODSEP), 

Propagating P by integrating P is more mathematically 
accurate than the use of effective process noise as in the J method, 
and it lends itself to greater modeling flexibility for Q* The draw- 
back to this method is the increased run time* 
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5.0 TRAJECTORY GENERATION - TOPSEP 

The targeting and optimization mode, TOPSEP, generates a reference 
trajectory which is supplied as basic input to the error analysis and 
simulation modes. The primary purpose of TOPSEP is to incorporate in 
this trajectory all of the desired flight characteristics for a particular 
interplanetary or near-Earth mission while optimizing the final spacecraft 
mass. Injection conditions (See Appendix 5), a thrusting time history, 
and other control parameters are found which accomplish this optimization 
and yet lead to the required target conditions. The target constraints may 
be the final spacecraft state (cartesian or B-plane coordinates), final 
orbital elements, radius of closest approach, or other mission specifications 
which are listed in Table 5-1 and in the input section of the Users Manual . 


Control Parameters 

Target Parameters 

Performance 

Parameter 

Initial State and Mass 

Impact (B) Plane 

Final Mass 
(Payload) 

Thrust Magnitude 

Sphere-of-Inf luence Time 


Thrust Direction 

Closest Approach Conditions 


Thrust Times 

(Radius, Inclination, Time) 


Base Power Level 

Target Centered State 


Exhaust Velocity 

Heliocentric State 



Table 5-1. Control, Target, and Performance Parameters 


The manipulation of trajectories to satisfy mission requirements is 
managed in three submodes of TOPSEP which represent successive stages of 
trajectory development. These submodes are: 
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1. Nominal trajectory propagation 

2. Grid generation > 

3. Targeting and Optimization 

a. Trajectory targeting only 

b. A combination of trajectory targeting and optimization 

c. Trajectory optimization only 

Generally, these submodes are employed in order as listed above. However, 
any submode may be skipped or used individually if the proper control pro- 
file is available. Due to the simplicity of the first two submodes a 
brief discussion of their operational procedures is all that is necessary to 
understand their analytical basis in TOPSEP. The targeting and optimiza- 
tion submode will be reviewed in greater depth. 

5.1 NOMINAL TRAJECTORY FWACATION 

The simplest TOPSEP application is propagation of a single trajectory 
for spacecraft ephemeris information. After all the trajectory parameters 
are initialized, the trajectory is propagated* from the initial state to 
the termination condition. TOPSEP performs no additional analysis of 
the trajectory when operating in this submode. This submode is also 
used for manual manipulation of the control profile. 

5.2 GRID GENERATION 

The grid generation submode is available to produce a number of tra- 
jectories which do not necessarily satisfy mission requirements but pro- 
vide a range of trajectory solutions. Thus, the main purpose of the grid 
submode is to locate desirable control regions for further examination. 


t 

/ 
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In turn, each control is incremented a fixed amount while the remaining 
controls maintain their nominal values, A single low thrust trajectory 
is generated for each control change and the associated target error index 
is calculated. Subsequently, contours of constant target error may te 
plotted in the control space so that some control regions can be eliminated 
from further consideration. Upon completion of the grid, the trajectory 
generation mode is terminated and the program user must choose the best 
control profile to initialize targeting and optimization or to employ 
another grid approach. 

5,3 TARGETING AND OPTIMIZATION 

The trajectory targeting and optimization submode features a discrete 
parameter iteration algorithm which accommodates the non-linear aspects of 
the low thrust problem. The algorithm is a modification of Rosen’s pro- 
jected gradient method (PGM) for non-linear programming (Refs. 5 and 6). 

The t ^rameters which have been chosen to shape the trajectory (Table 5-1) 
constitute the control profile and sre subject to modification by the PGM 
algorithm. Based upon the sensitivities of the final S/C mass and state 

to control variations, corrections to the profile are computed which 

» 

maximize performance while minimizing target errors. The performance is 
measured simply by the value of the final spacecraft mass while the target 
errors are measured according to the constraint violations. The method 
chosen to represent the target errors in terms of a scalar measure is the 
quadratic error index which is the weighted sum of the squares of the 
target errors. 

When the targeting and optimization submode is entered, a nominal 
trajectory is propagated directly from the input parameters. A series of 
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tests is performed to determine which submode (targeting, optimization or 
both) is to be executed. If the target error index is large, the submode 
will be exclusively targeting. However, a target error index smaller than 
some value (TUP in namelist $T0PSEP) will result in simultaneous targeting 
and optimization. Whenever the index is below a specified lower bound 
(TL0W in namelist $T0PSEP), the optimization algorithm will be executed. 

Prior to the application of the projected gradient algorithm, the 
targeting sensitivity matrix S and the performance gradient £, are computed. 
Elements of the S matrix represent the sensitivities of individual target 
parameters to changes in controls and are used for both targeting and 
optimization. Similarly, the elements of the £ vector represent the 
sensitivity of the performance index to changes in controls although 
these elements are used only for optimization. For purposes of targeting 
only, S is computed from the integrated state transition matrix (STM) 
and £ is ignored. Appendix 7 discusses the formulation of S from the 
integrated STM. Whenever optimization is to occur both S and £ are con- 
structed by finite differencing techniques. Following the determination 
of S and £ a weighting matrix which amplifies or diminishes the effects 
of the chosen controls is calculated. Applying the projected gradient 
algorithm a control correction is established. The magnitude of the con- 
trol change is determined by computing trial trajectories. The new control 
profile is simply the old control profile plus a scalar multiole of the 
control correction such that the targeting error index is minimized and/or 
the performance index is maximized. If the optimization is complete (the 
values of the performance index have converged to a maximum), TOPSEP is 
terminated. Otherwise, the submode decision is made again and the cycle 


is repeated 
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Of primary importance in the targeting and optimization submode is 
the selection of the control correction. In the following sections this 
selection process will be discussed, 

THE PROJECTED GRADIENT METHOD 


The projected gradient method has been devised to maximize a per- 
formance index while simultaneously minimizing an error index. Since 
maximizing performance is equivalent to minimizing cost in an optimization 
sense, PGM's purpose relative to the trajectory problem may be restated as 
minimizing fuel expended as well as minimizing target error. The concept 
of net cost, which is simply a more realistic assessment of fuel expended, 
will be discussed later in this section. 

The projected gradient algorithm employs cost-function and constraint 
gradient information to replace the multi-dimensional targeting and opti- 
mization problem by an equivalent sequence of one-dimensiona 1 searches 
(Ref. 7). In this manner, it solves a difficult multi-dimensional problem 
by solving a sequence of simpler problems. In general, at the initiation 
of the iteration sequence, PGM primarily satisfies the constraint require- 
ments. As the iteration process proceeds, the emphasis changes fron con- 
straint satisfaction to cost-function reduction. 

Since numerous analytical developments of this technique are available 
(Refs. 5 and 6), this presentation will primarily emphasize the geome- 
trical aspects of the algorithm. Clearly, the geometric interpretation of 
the algorithm is the motivation for the logic contained in TOPSEP, and a 
basic understanding of these concepts is usually sufficient to enable the 
user to efficiently manipulate ^TOPSEP input and to handle diverse mission 
problems. 


i 
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PROBLEM FORMULATION 

The projected gradient method solves the following non-linear pro- 
gramming problem: 

Determine the values of the control variables, u 9 that minimize the cost 
function (optimization variable) 


F(u) 

subject to the equality constraints 

e (x(u)) = T (x(u)) - l d = 0, 

u £ U, an M-dimensional control space 
T £ T, an N-dimensiona 1 target space 
x £ X, a six-dimensional state space 
M > N 

where the elements of £, T , T^, and x are referred to as the target error, 
the target values, the desired target values, and the final state respec- 
tively; and F is a scalar valued function measuring system cost. 

In an attempt to solve the constrained optimization problem, iterative 
methods are employed. The following scheme briefly describes the process 
occurring in TOPSEP. 

• Guess u (In general u will not satisfy the constraints nor 
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AO 


• Determine a correction fan such that 

•• F(u o + Au) ^ F (Uq) and 

•• || §.( +AE)|| < II e <1^)11 

• Iterate until 

•• F is minimized and 

•• II ill < C, a pre-determined tolerance 
NOMENCLATURE AND CONCEPTS 

To fascilitate the discussion of the projected gradient algorithm, 

the following nomenclature and basic concepts will be introduced, x denotes 

• column vector whose elements are x^, where i = 1,2, •••, n and n is the 

T 

dimension of the space containing x. Y denotes the transpose of the real 

matrix Y- The feasible region defined in the M-dimensional control space 

within which PGM operates is the restricted space 

V i u, A u , i ■ 1, 2, ***, M. 

N 1 MAX 

The equality condition implies that the control is on a bound. The cost 
gradient g is an M-vector of partial derivatives and is defined as 



The sensitivity matrix is that matrix whose rows are the gradients to the 
equality constraints, and is denoted by 

S(») - 9 1 , 

a a. 

where £ is an N-dimensional vector. The target error function is defined 

to be 

E(u) - e T [tT] e 



where is a target weighting matrix which will be defined later. 

Corresponding to each control vector £ in the control space U, 
there is an error vector e. Let A be the set of all u such that 

"" O “ 

e (u) ■ 0. 

A then represents all the control vectors satisfying zero-target error, 
o 

It can be shown (Reference 6) that A^ defines an M-N dimensional non-linear 

hypersurfcce or manifold in U. Unfortunately, A^ cannot be defined 

explicitly; hence, one cannot easily find a u which is an element of A q . 

However, A^ can be estimated implicit'/ via the sensitivity matrix. 

Let A be the set of all u such that 
c ~ 

e(n) ■ 

where c is a vector of constants. Thus, A represents a non-linear mani- 
— c 

fold containing those control vectors which provide constant target error. 

It can also be shown (Reference 6) that any u in the control space is 

contained in one and only one A^, At a given u, the corresponding 

non-linear manifold A may be approximated by a linear manifold B(u) which 

c 

is defined explicitly by the sensitivity matrix S(u). The linear manifold 

B(u) may be considered a tangent hyperplane to A at u. The orientation 

c 

of B(u) in the control space allows one to define a search direction to A 

o 

which is orthogonal to B(u). This search is in the direction of maximum 
decreasing target error. 

Let B(u) denote the orthogonal complement to B(u). One can demon- 
•irate (Reference 6) that B(u) la the unique linear space that can be 
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translated to obtain the linear manifold B(u). Furthermore! there 
exist unique orthogonal projection operators P(u) and P(uJ that resolve 
any vector In the control space Into Its corresponding components In 

ft# 

B(u) and B(u), respectively; that is 

A# 

u ■ P(u) u + P(u) u. 

In particular, 

P(u) “ S T (SS T ) S (5-1) 

and 

P(u) - I -*P 

Al 

where I is an identity matrix. The projection operators p and P thus pro- 
vide a simple method for reconstructing a general vector Au emanating from 

A# 

u into its two components in B and B. In a discussion to follow later in 

e* 

this section, it will be explained how Au^ may be defined such that PAo 
represents he .ertrol correction to minimize the target error and PAtj 
represents the control correction to minimize cost. 

The final key concept employed by PGM is the idea of problem scaling. 
The purpose of problem scaling is to increase the efficiency of the tar- 
geting and optimization algorithms by transforming the original problem 
Into an equivalent problem that is numerically easier to solve. 

To numerically scale a problem, two general types of scaling are 
required: 1) control variable scaling, and 2) target variable scaling. 

Control variable scaling is accomplished by defining a positive diagonal 
scaling matrix, (UWATE in namelist $T0PSEPX such that the weighted 
control variables are given by 

u* - W u . 

— u •“ 
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Similarly, target variable weighting is accomplished by defining a 

positive diagonal scaling matrix, W (TART0L in namelist $T0PSEP) , such 

e 

that the weighted target variables are 

T' (u) * fw *] T (W u') . 

— — lei" u — 

The target error index is then 


E(u) 


„,T 
e e 


TOPSEP contains several options for computing the control variable 
weighting matrix. The option most often used is the normalization scaling 
matrix (See Appendix 6 for other options). 


hi 


ii 




The target variable weighting matrix is always computed as the reciprocal 
of the constraint tolerances and is given by 

hlu- -}-• 


, th 


where is the tolerance for the i target error. 

For simplicity, the following discussion of the algorithm assumes an 
appropriately scaled problem. However, the scaled equations can be obtained 
by making the following simple substitutions. 

u replaced by u/ 
e replaced by e 1 


S replaced by MOlhl' 1 
£ replaced hi ' 1 
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DIRECTION OF SEARCH 

The concept of the direction of search in control space needs 
•lightly more elaboration. The direction of search is nothing more than 
• particular line in the control space along which the target error is 
reduced or along which the cost function is decrease!. In a more precise 
tense, the direction of search at i is a half-ray emanating from u^ Thus, 
for any positive scalar, V , the equation 

u “ (j +)£au, 

sets the limits of this half-ray and represents a "step" in the direction 
Au from u^ . This is illustrated in Figure 5-1. 



Figure 5-1. Direction of Search 


This concept of direction-of-search is particularly important since it 
enables the M-dimensional non-linear programming problem to be replaced by 
a sequence of one-dimensional minimizations. What remains to be explained 
la: 1) how to select the direction of search, and 2) how to determine 

the step size in that direction. 

i 

u 


1 



1 


- 1 
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The projected gradient method u®es two basic search directions. 

For this discussion they will be termed constraint and optimization 
directions. PGM proceeds by taking successive steps in one or the other 
of these two directions. The computation of each of these search direc- 
tions is described below at a particular point u in the M-dimensional con- 
trol space where N constraints (target conditions) are enforced. 

CONSTRAINT DIRECTION 

Hie constraint direction depends critically on the number of targets. 
Two cases are distinguished below: 

I. If N<M, then that unique control correction 4u is sought which 
solves the linearized constraint equation 

^S^Au + e (u) ■ 0 (5-2) 

and minimizes the norm of » The solutions to the preceding 
vector equation define the M-N dimensional linear manifold B q (u) 
which is an estimate of the non-lioear manifold A q (zero-target 

A 

error). The desired minimum norm correction A u is then the vector 
of minimum length in the control space from vi to the linear mani- 
fold B (u), thus requiring Au to be orthogonal to B (u). 
j ~ o 

Application of the linear operators P and P allow one to represent 

db ^ 

Au as the sum of two orthogonal vectors relative to B (u) or 


A *** • A 

Au ■ P Au + P A u ; 


however , 




P Au - 0 

since there are no components of Au in B^. From equation (5-1) 
A u. may be expressed as 


r 

1 



which may be reduced to 

Au - -S T (SS T ) e(u) ' 

using the constraint equation (5-2). This correction is illus- 
trated in Figure 5-2. 


£u f minimum norm 
correction 


B (u), intersection of 
o *■" 


linearized constraints 



Figure 5-2. Illustration cf Minimum Norm Constraint Direction 
for N - 2, M * 3. 

The direction of search then is simply taken to be this minimum 
norm correction to the linearized constraints. 

If N ■ M there is unique solution to the linearized constraint 
equations without the additonal requirement that the norm of 
the control correction be minimized. The solution for Au_ reduces 
to the familiar Newton-Raphson formula for solving M equations 


with M unknowns; namely 



I 


- s * e (u ) 


The Newton-Raphson correction is illustrated geometrically in 


Figure 5-3. 


Third linearized 
constraint (e^ = 0) 


/m" 










First linearized 
constraint (e^=C) 


life 


Second linearized 
constraint (e^O) 


Ay. > Newton-Raphson 
correction 


W , ' /. ■ V 






[■ / / 7 '///.%' 


B (u), intersection of 
o 

linearized constraints 


Figure 5-3. Illustration of Newton-Raphson constraint direction for 
N « M - 3 

OPTIMIZATION DIRECTION 

When the number of targets is less than the number of controls, it 
is then possible to minimize the cost function F (u ) assuming, of course, 
that li is some non-optimal control profile. Obviously, the steepest 
descent direction, - g(u), would be the best local search direction for 
reducing the cost function. Such a direction, however, could produce 
unacceptable constraint violations. To avoid this difficulty PGM ortho- 


gonally projects the unconstrained negative gradient,- g, onto the local 
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linearized constraint manifold B (u) . By searching in the direction of 

c 

this negative projected gradient the algorithm can guarantee in a linear 
sense that there is no further constraint violation than tha'~ of e,(u) , 

To calculate this direction, it is only necessary to apply to - g the 
projection operator P(u) which will map the vector into its component on 

aw 

the linear manifold B c (u) . Thus, 

A A 

Au = - P g (u) 

- -t 1 - p]g (h) i 

- -[i - S T (SS T ) s] g (S) 

COMBINED TARGETING AND OPTIMIZATION DIRECTION 

When it is desirable to minimize the cost function as well as reducing 
the target error the constraint direction and optimization direction may 
be combined such that the resulting control correction is of the form 

Au ■ AUj + A 

where Au^ is the optimization correction and Au^ is the constraint correc- 
tion. Note that andAu^ are orthogonal components of Au . Depending 

upon the magnitude of the target error, one may want to emphasize either 
optimization or targeting. Since this decision is rather subjective and 
linked directly to the degree of problem non-linearity, an option is pro- 
vided to weight Au (See Page SO) in one of its component directions. 

Figure 5-4 and Figure 5-5 illustrate the geometric interpretation of the 
resulting control correction. 


\ / 


Nonlinear Constraint 
Manifold (e=0) 


Constrained 

Optimum 


N 


Figure 5-4. Geometric Interpretation of a Combined Targeting 
and Optimization Control Correction, N=l, M=3 



B c (i) 


Figure 5-5. Illustration of Combined Targeting and Optimization 

Control Correction As Seen In *B (u), 

c 

y A 

(Enlargement of B^(u) from Figure 5-4). 
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The total control correction is constructed as follows, where d is an input 
scalar. 


A u 


1 

UuJ 

1 

IK 1 

1 


* d * AUj + Au^ 


(DP2 in namelist $T0PSEP) . Thus , one has the flexibility of determining 
the magnitude of the optimization correction relative to the magnitude of 
the corstraint correction. The optimization correctiun can then be written 


as 


= -Ji (SS T )" e * (1 + d 2 ) 


(I-S T (SS T ) S) g . 


|| (I - S T (SS T ) ^s, g || 


The norm of the control correction Au^which is obtained by sumraingAu^ and 
/i not as important as the direction. The resulting half-ray provides 

the basic search direction with which to calculate the trial step. 

TRIAL STEP-SIZE CALCULATION 

At any particular point u in the control space, the PGM algorithm 
proceeds by reducing the multi-dimensional problem to a one-dimensional 
search in the direction prescribed by the constraint or optimization con- 

A 

trol change vector. Once the initial point u and the direction of search 

A 

&u are specified, the problem reduces to the numerical minimization of a 
function of a single variable, namely , the step scale factor V . PGM per- 
forms this numerical minimization by polynomial interpolation based on 
function values along the search ray and the function’s value and slope at 
the starting point. 

CONSTRAINT DIRECTION 

A 

The function to be minimized along the constraint direction Au^ is 
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E($), the sum of Che squares of the target errors. 

E(fc) " £. T (u+Au 2 ) « (u + a H2> 

Evaluation of the function at l » 0 results in 



Differentiation via the chain rule yields 


a# 


•J* A 

= 2 e (a) s A 


¥= 0 


* 



If constraints are reasonably linear, a good initial estimate for the 
minimizing & is one. 


OPTIMIZATION DIRECTION* 


The function to be minimized along the optimization direction & u^ is 
the estimated net cost function C(Y), where 


C(» 


) «= F(C +*4^) - F(u) + g - S T (SS T ) 1 e (f+U^)*] 


Change in cost produced 
by a step of length 
% || A Hl || along 


Linearized approximation to change 
in cost required to perform a minimum 
norm correction in order to retarget 


Clearly, 

C(0) = - g sT (SS T ) e <£) 

By expanding C(0) in a taylor series in 't about * 0, and by making use 

A 

of the fact that PAu^ * 0 , it can be shown that 


3c(* ) 


d* 


* UE, 


TS * 0 

These properties are illustrated in Figure 5*6. 


Both the constraint and optimization directions are based on a 
sensitivity matrix assuming a linear space. Due to nonlinearities in the 
problem, it is often necessary to restrict the trial step such that 
unexpected increases in cost O’- target error are reduced. In addition, 
the control correction Au does not reflect the "nearness" of control bounds 
as long as u. is not on any bound. Thus, the trial step must also be 
restricted so that the new control vector remains within the bounded con- 
trol space. For these reasons a maximum limit is placed on IS . After 



Figure 5-6 


Properties of net-cost function along the direction of 
search. 
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Au is computed, a calculation is made to find the 'C .along the 

\ 

search direction which will size a trial step that: intersects the 
boundary of the feasible region. This value is compared to the 
maximum step allowed by the user to counter the nonlinearity prob- 
lem. The smaller value is specified as the maximum IT allowed 
during the search. A method has been devised to alleviate the 
problem of a control which is very near a boundary. A tolerance 
region is defined in the neighborhood of the boundary surface such 
that if the control is within this region and Ah intersects the 
boundary the search to minimize the net cost or target error can 
continue along the boundary without calculating a new sensitivity 
matrix. Once a control element reaches one of its bounds it becomes 
inactive. Unless a subsequent correction for this control element 
is back into the feasible region it remains inactive. 

ONE DIMENSIONAL MINIMIZATION 

Nonvariant minimization in PCM is performed exclusively by 
polynominal interpolation. The actual function to be minimized is 
fitted with one or more quadratic or cubic polynominals until a 
sufficiently .'ccurate curve fit is obtained. The minimum of this 
curve and the corresponding scale factor can easily be found analyt- 
ically. 

The one-dimensional search proceeds by taking trial steps in 
the Ah direction to obtain information about the function to be 
minimized. If t*Au is a constraint correction, the quadratic 
error function is evaluated; if TfA H is an optimization correction. 


i 
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the net-cost function ij evaluated; and if TTAij is a combined cor- 
rection^ function which is a weighted combination of the error 
function and net-cost function is evaluated. 

The minimization routine makes ingenious use of all the infor- 
mation it accumulates. The following curve-fitting techniques are 
applied in order. 

1. Quadratic polynomial fit: two points-one slope; 

2. Cubic polynomial fit: three points-one slope; 

3. Quadratic polynomial fit: three points; 

4. Cubic polynomial fit: four points. 

Each time a trial step is taken, the function whic.i is evaluated is 
used as a trial point to analytically determine the next trial step. 
The analytical formulation of the above curve fits may be found in 
the subroutine description of M1NMUM in the Program Manual . 


! 

i 
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6. LINEAR ERROR ANALYSIS - GODSEP 


GODSEP analyzes spacecraft and trajectory related dispersions as 
a function of expected uncertainties in dynamic and navigation parameters. 
The ensemble of expected errors is studied without actually simulating 
individual trajectories by applying linear techniques. That is, small 
deviations about a reference trajectory are linearly related to other 
deviations by a transformation matrix. For example, the state transition 
matrix relates position and velocity deviations about the reference tra- 
jectory from one time point to another. The ensemble of errors, or covar- 
iance, is assumed to have a zero-mean Gaussian distribution, except for 
special processes. 

Probabalistic a-priori errors in the environment, spacecraft and 
tracking systems are propagated in time along the reference trajectory 
through sequential events such as orbit determination (OD) and guidance 
corrections. Two types of ensemble error or covariances are distinguished 
knowledge , which reflects the ability of the OD algorithm to estimate the 
spacecraft state and other parameters; and control . which represents the 
dispersions of the actual spacecraft trajectory about the reference. 
Covariance propagation is done by either integration of covariance varia- 
tional equations, or by the state transition matrix method. 

The error analysis proceeds sequentially from start time through each 
specified trajectory event to final time. Event types available are 
measurement, propagation, eigenvector, prediction, thrust switching, and 
guidance. A measurement event processes tracking data at a time point by 
applying the user specified OD algorithm. Available to the user are both 
Kalman-Schmidt (K-S) and sequential weighted least squares (WLS) filters. 
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GODSEP modularity a^jo allows the user to insert his own filter algorithm 
quite easily. The filters are distinguished by their methods of gain 
matrix calculation and subsequent update of the knowledge covariance. 

A propagation event merely updates the knowledge (and control) covari- 
ance at the event time. Its primary value is in maintaining accurate covari- 
ance values during long propagations by forcing computation of the effective 
process noise over predetermined, user-specified intervals. 

An eigenvector event is used for information display and behaves similar 
to a propagation event. Covariance matrix sub-blocks are converted to 
standard deviations and correlation coefficients. It also computes eigen- 
values, their square roots, and eigenvectors for the position and velocity 
3x3 sub-blocks of the state covariance matrix. Thrust switching events are 
simply eigenvector events at the time where a change in the number of 
thrusters or thrust policy has occurred. 

A guidance event is an update of the control covariance to reflect 
Implementation of a trajectory correction. A correction is not performed 
deterministically, but only in a probablistic sense. The guidance event 
computes expected correction covariances (Av or thrust control), target 
error covariances before and after the guidance event, and the updated state 
control covariance. 

The following sections will describe in more detail the analytical 
foundations of GODSEP. 

6.1 AUGMENTED STATE 

The augmented state discussed previously in TRAJ (Section 4,5) 

Includes dynamic parameters besides the basic spacecraft position and 
velocity vectors. In GODSEP, the augmentation not only adds measurement 
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related parameters to this list, but also distinguishes between solve- 
for and consider. Solve-for parameters arc directly estimated by the OD 
process. Consider parameters are system uncertainties which are recog- 
nized and accounted for in the estimation algorithm but are not estimated, 
usually because the process cannot be adequately modeled or there is a high 
correlation between two (or more) parameters which might cause numerical 
difficulties if both were solved-for. 

The possible augmented parameters that can be either solved-for or 
considered are 

thrust bias (magnitude and pointing) 

position and velocity of a selected planetary (ephemeris) body 
gravitational constants of ephemeris body and/or sun 
tracking station locations 
sensor bias (range, range-rate, etc.) 

Time varying thrust noise (magnitude and pointing) can only be considered 
in the standard CODSEP analysis, but can be -solved-for (or considered) in 
the covariance integration option (PDOT in Section 6.2). A third possible 
category, in addition to solve-for and consider, is the ignore parameter 
used in generalized covariance analysis (Section 6.5). 

The total ensemble of state uncertainties, or error covariance, 
Including all augmented parameters, is formed by applying the expected 
value operation on state deviations from their reference values. 

»-t[fcSs T ] 

, The covariance P contains uncertainties and their respective correlations 

for all parameters in the augmented state. There are two covariances, 

U 


dynamic 




measure- 
ment 


{: 
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KJ 


corresponding to contxol and knowledge, which are computed in parallel 
through a GODSEP analysis. Starting with a-priori values, each P is 
modified between events by trajectory propagation effects, and at events 
by either OD (knowledge) or by guidance corrections (control). 

6.2 COVARIANCE PROPAGATION 

There are two methods available in GODSEP for propagating covariances 
between events: transition matrices ($) and explicit covariance integration 

(PDOT). Although these two techniques were discussed in TRAJ, Sections 4.5 
and 4.6, respectively, their importance in GODSEP requires additional explana- 
tion. 

The most coimon form of covariance propagation, both in GODSEP and in 
other linear error analyses, makes use of transition matrices. This is 
because the 4' s are a characteristic of the trajectory, not of the covari- 
ance, A covariance P(t^) is propagated from time t^ to, t^ by 

• P(t 2 ) - rf 21 PC^) <t 2 \ + Q 21 (6-1) 

where ■ 4 ( t 2 » t 1 ) and Q 21 “ Q(t 2 ,t^) is an effective process noise 
covariance. 

Transition matrices can be stored, for example on magnetic tape, to 
be used for analyses of different error source levels and navigation strate- 
gies. Obviously, if ^'s have been computed between the intervals t Q , t^, 

t«, *»*t , they can be used to propagate any P as long as the set of pro- 
2 N 

i 

pagation times is a subset or equal to the original set. Transition matrices 
can always be chained to cover desired propagation intervals, for example, 

letting ^ . ^(t 3 ,t x ) - 4(t 2 , t x ) 

then P(t 3 ) - 4 n P(t t ) 0 31 + Q 31 
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The method of computing and storing 4' s (on the STM file) over a grid 
of time points Is used in GODSEP to facilitate parametric error analyses. 

Since covariance propagation accounts for uncertainties in all 
dynamic parameters which have been augmented to the basic spacecraft state, 
the transition matrix must have the same augmentation. In actual operation, 
TRAJ provides a transition matrix containing only dynamic variations which 
GODSEP must augment with appropriate rows and columns of zeros (for measure* 
sent parameters) such that the total augmented 4 is consistent with P. 

An additional requirement for GODSEP is the modification of the thrust 
sensitivity matrix 9 computed by TRAJ as part of the augmented transition 
matrix. 




^x(t2) 

Tu 


where u are constant thrust controls (proportionality, cone and clock) over 
the Interval (t^,t^). The 9 matrix is used in GODSEP to map thrust biases 
into spacecraft state uncertainties. However, GODSEP chrust biases refer 
Co a single thruster. If more thrusters are operating, and if each operating 
thruster is assumed to be independent of all others in terms of bias, then 
Che total effective bias is simply the single thruster bias divided by 
where N is the number of thrusters, or 


•i - F % 

The effective process noise, Q , is a very important conditioning 
term on the covariance propagation. Because a rigorous mathematical com- 
putation of Q involves (1) modeling of a process or processes which are 
ill-defined and (2) evaluation of complex double Integrals, GODSEP uses a 


i 



f 
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simple analytic approximation. The effective process noise assumes that 
time-varying thrust errors appear as stationary first-order Gauss Markov 
processes. The moie rigorous modeling is performed in the PDOT option to 
be discussed shortly. The relationship between P and Q precludes the aug- 
mented state from containing time-varying thrust terms so that process noise 
takes on the app'-r'rance of consider parameters. 

The effective noise over a time interval t^ to t^ directly affects 
only the spacecraft position and velocity uncertainties at 


< 21 


1/2 At 


[C y 


+ l 


r°)i! 

•21 to *21 


where J is the 6x6 state transition sub-matrix of the augmented tranrition 
matrix (4), 

H, 


g(tj> 


8 (t t ) 


H. 


g(t 2 ) 0( g (t 2 ) 
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1 


variances in thrust proportionality, cone, clock, 
respectively. 

nu. • er of operating thrusters 

(see Section 4.1 and Appendix 4) 

a** 


T l* V 1 3 


process noise correlation time# 
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Thrust proportionality error is scaled by the number of thrusters 
because it is assumed that time-varying noise is independent for each 
thruster, just as bias is. Thrust pointing noise is also scaled 
because it is assumed to be caused by the thrust vector control system 
which currently consists of gimbaling each thruster independently. 

This empirical model for is generally effective over propagation 
intervals on the order of 50 days or less. Propagation events can be 
employed in GODSEP to break up longer intervals and to ensure the accuracy 
of Q. 

The second method of covariance propagation, PDOT, is used primarily 
to examine thrust noise effects. Process modeling is mathematically 
rigorous and includes augmentation of thrust noise parameters to the basic 
state. Recalling from Section 4.1, the linearized equations of motion, 

£ x = F S x 

and from Section 4,6, the corresponding covariance matrix differential 
equations 

P * FP + PF T + Q 

where x is the augmented state which, for PDOT, includes at most space- 
craft position and velocity, thrust biases, thrust noise and tracking sta- 
tion locations, F is the variation matrix, and Q is a white noise term 
which affects only the thrust noise directly. If thrust noise is omitted, 
then the integrated covariance would ir. theory be identical to a similarly 
augmented covariance propagated by transition matrices. 

In PDOT, the time-varying noise is modeled as a stationary Gauss- 



Markov process, as in Q, 


r 


0 


cJ = 


- 1 
T, 


6 4 


+ W 


where is a 6 x 1 vector of independent noise parameters corresponding to 

thrust proportionality, cone and clock, each of which is described by two 

processes having their own distinct correlation times (T) . This permits 

the study of superimposed and multi-process effects. W is a white noise 

component which drives the time varying noise (and defines the only non- 

T 

zero term of Q which is E [w W ] ). 

Since (O is in the desired form of the linearized equations of motion, 
it can easily be augmented to the state vector (and covariance). Thus, 

60 can be solved-for in the PDOT mode although in reality this practice is 
questionable because of the w modeling assumptions - who knows how thrust 
noise really behaves? 

One of the more useful applications of PDOT is in refining the form 
of effective noise, Q, for a particular mission and in verifying the 
explicit assumption in Q of zero correlation between noise and state para- 
meters . 

Whether state transition matrices or PDOT is used for covariance 
propagation, an auxilliary computation is the vehicle muss uncertainty. 
Since mass and thrust magnitude uncertainties ate indistinguishable 
ir their trajectory effects, that is, they are correlated one to one, 
GODSEP has chosen to model thrust (accelerat ion proportionality) magni- 
tude explicitly, and provide the approximate equivalent mass uncertainty 
as supplementary information. 
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Two types of mass uncertainty are distinguished: knowledge 

(estimated) and control (actual). Estimated mass uncertainty is the 
Instantaneous knowledge error in thrust magnitude, bias and noise. 


estimated (J* ^ * (CTu^ + (T*^)* r ‘ 

m ab an 

2 2 2 

where (T , (T., (T are the variances in mass, thrust bias 

m ’ w ab w an * 

proportionality (from P^) and thrust noise proportionality (from 
or Q), respectively. 

Actual mass uncertainty is the cumulative mass variation reflected 
by the control error covariance. The actual mass deviation from the 
reference at time t + A t based upon uncertainties from time t is 


actual <rj (t + A t) * f ff^Ct) + * 





+ 2 «*.„ 2 t 2 * 


where m, and CT fln are averaged ever the interval At, c is the 

exhaust velocity and X the correlation time. <J^ b (and CT r for PDOT) 

are obtained from the augmented control covariance. Accuracy of the 

mass variance computation depends upon the event schedule because GODSEP 
2 

evaluates (JT from one event to the next, 
m 
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6,3 MEASUREMENT TYPES 

In a linear error analysis, the reference trajectory deterministically 
characterizes the motion of the S/C, and no real state vector estimation 
is explicitly performed in GODSEP. Rather, an orbit determination anal- 
ysis estimates how well the state vector can be determined if the S/C 
were to move along the reference trajectory and were to be observed as 
directed by the analyst. In this sense, the term "orbit determination" 
refers to the calculation of a knowledge covariance based upon the proc- 
essing of modeled observational data. This section of the Analytic Manual 
describes the data types and mathematical models that have been implemented 
in GODSEP. The next section will treat the problem of filter formulation 
and the process of updating the knowledge covariance. 

When an observation is to be simulated in GODSEP, the knowledge 
covariance is propagated to the scheduled measurement point and is made 
available to be updated by the filter. Before this can happen, it is 
necessary to evaluate the observation matrix which relates the observables 
to the state vector. Given an arbitrary veci - (or scalar) measurement 
Z ~ Z (X) w here X is the total augmented state insisting of 


► « 

X 


spacecraft position and velocity 

£ 


solve-for parameters 

U 

s 

dynamic consider parameters 

V 


measurement consider parameters 

w 

* ““ * 


ignore parameters 


then the linearized measurement, which assumes small deviations from the 
nominal, is 

S i * 


H $x+H Ss+H Su+H & V r H S S 
X u - v w 
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where H 


2 .'- H 

9 ) • • i | n 

X w 


^ z 


are the observation matrices, 


X X ' ' w ^ w 

all of which are computed analytically in GODSEP. 

GODSEP has the capability of processing the following measurement 


types : 


in- 


earth- 

based 




spacecraft 

based 


{ 


o 2-way range 
o 2-way doppler 
o 3-way range 
o 3-way doppler 

o simultaneous 2-way and 3-way range 
o simultaneous 2-way and 3-way doppler 
o differenced 2-way and 3-way range 
o differenced 2-way and 3-way doppler 
o azimuth and elevation angles 

o right ascension and declination angles (of target 
body) 

o star-planet/target body angles 
o planet limb angles (apparent planet diameter) 


All earthbased data types which nmke observations of the S/C are 
applicable to both near earth and deep space missions; however as a 
practical matter, azimuth and elevation angles are normally used for 
near earth analysis only. Astronomical observations of the apparent 
right ascension and declination of the target body are used to determine 
ephemeris errors and can be made concurrently with earthbased tracking 
of the S/C. (Astronomical observations provide information about the 
state of the target body and indirectly imply information about Lne S/C 
motion if there are dynamic and/or measurement correlations between the 
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S/C and the target body). Spacecraft based observations have been 
formulated to model measurements of the target body with an onboard 
optical system. Hence, GODSEP permits an integrated navigational 
analysis for interplanetary missions with astronomical, radio, and 
onboard optical measurements all in one computer run. 

For earthbased radio observations of the S/C, the program normally 
uses the standard DSN tracking stations located at Goldstone, Madrid, 
and Canberra. However, the locations of these stations may be changed 
by input and as many as six others added. For astronomical observations, 
the nominal observatory is assumed to be Kitt Peak in Arizona. As with 
the DSN stations, the location of the observatory may also be charged to 
model whatever real observatory the analyst chooses. 

In order to display the analytic partials contained in the observa- 
tion matrices, earthbased data types are separated into two categories 
according to their normal application for deep space or near earth 
missions. Thus, all mathematical models for range and range-rate data 
will be described as deep space data, although they are directly appli- 
cable, without reformulation, to the near earth problem. The right ascen- 
sion and ^lination observations are also grouped in the deep space cate- 
gory although visual observations of near earth satellites could be made 
with a few minor changes to the model. Only azimuth and elevation angle 
measurements are specifically treated as near earth data, and a discussion 
of their observation paitials will follow the development of the deep 
space data types. Finally, a description is given for the S/C based data 
types and the observation partials. 



1 


64-C 

The following definitions of position and velocity vectors are neces- 
sary to relate the mission geometry to the observable quantities. All 
vectors are assumed to be column vectors and are expressed relative to an 
inertial, ecliptic coordinate frame, unless otherwise noted. 


x, x 


S/C heliocentric cartesian position and velocity 



Earth heliocentric cartesian position and velocity 




Target body heliocentric cartesian position and 
velocity 


— i*— i 


Station 1 geocentric cartesian position and velocity 


>*2 * Station 2 geocentric cartesian position and velocity 

- S/C position and velocity relative to Station 1 

» 

P_2’^2 = S/C position and velocity relative to Station 2 


Unit vectors defining direction of S/C from 
Stations 1 and 2 respectively 


* * 

P l - P 2 





'-2 


S/C range and range-rate from Station 1 
S/C range and range-rate from Station 2 
3-way range and range-rate 

Differenced 2-wav and 3-way range and range-rate 

Geocentric cylindrical coordinates of Stations 

1 and 2, s * (r , , r) T 

s 

Zero vector, 3x1 

Identity matrix, 3x3 unless noted otherwise 



Figure 6-1. Tracking Geometry for Range and Range-Rate 
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We first note one identity which is used numerous times in the follow- 
ing derivations. Given the vector a. * b. * £ and its corresponding 

unit vector 3 • a/ |aj , 




ITT 1 1 


* * T ^ 


a a 


( 6 - 2 . 1 ) 




Two-way range and range rate from Station 1 are modeled 


IT ^ 

fi ■ $ <°l 


where 


n ■ s p i 




P ... 

-1 - ■ *g - % 


Differentiation yields the following results. 


3p yw - e. T ap / a , + V *%*. 


- £l • -5- [ 1 - Pi V ] + K • 
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( 6 - 2 . 2 ) 


The remainder of Che partials are produced in like manner but for 
brevity only the results will be printed 







£(*! *}) 



5 (Sj » ) 

3 % 


( 6 -?. 3 ) 



(6-2.4) 


(6-2.5) 


For use in Equations 6-2.3 and 6-2.5 above, we need the dexivative of 
the instantaneous geocentric ecliptic cartesian state of the tracking 

a* 

station with respect to its cylindrical equatorial coordinates of spin 


i 
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radius (r ), longitude (JO and 2 -height ( 2 ). If G represents the 

S 

Greenwich hour angle at launch, t the universal time (U.T.) at the 
launch epoch, t the U.T, at the current epoch, and 63 the Earth's 
sidereal rotation rate, we have 


*1 " E 


^ r g cos [ X + G + 6* (t-t o > ] 
r sin [ X + G + CO (t-t )] 


Xj « E 


- CJ r sin £ X + G +6J(t-t )] 

8 O 

6) r g cos [ X + G +40 (t-t Q )^ 


where E represents the 3x3 transformation from geocentric equato- 
rial to geocentric ecliptic cartesian coordinates. This transformation 
is assumed constant. Even though the Earth's obliquity to the ecliptic 
does vary slightly, its effect is negligible over the duration of the 
missions for which these programs are used. Letting 


X + G +» (t-t o ) 


0 



sin 0 


0 




, z) « E 


cos 0 -r 


sin 0 r cos 0 0 


( 6 - 2 . 6 ) 



-CJ sin 0 

*<ur cos 0 
8 

0 

CJ cos 0 

-Cj r sin 0 
s 

0 

0 

0 

0 


(6-2.7) 


Three-way range and range-rate are measured with one station on the 
DSN uplink and another station on the downlink. Three-way data may 
be processed by itself, simultaneously with conventional two-way data, 
or as differenced two-way minus three-way data, also known as QVLBI 
(quasi-very long baseline interferometry). Three-way data types are 
modeled as the sum of the two-way types plus a timing error term for 
ranging and a frequency bias term for range-rate. 

f 3 - + P 2 + c tt 

0 • p + p + c -Ax 

r 3 m *2 c f 

where Ac is the timing error, c the speed of light, and ^ / t the 
frequency bias term which results from drift error between the frequency 
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standards at the two separate tracking stations. The sensitivity 
partials for the three-way data types are formed by adding the par- 
tlals computed for each station individually. The c A t and 
c ^*/f terms arc treated either as biases or part of the white 
noise term. The differenced data types are modeled: 

Ap - - ? 2 * at 

A9 • K - P 2 - c if/£ 

The partials for the differenced data types are formed by differ- 
encing the individual partials, with the following exception. Since 

a a? r* t 

. irl. - 0 - P 

dx 3 X L v 1 v 2 

a w' 1 Vj and are very nearly equal (as are and £,,) ^ or 

i:terplanoc 2 t / missions, we use the following substitutions 
» 

£ x “ x _2 “ 

AC - Lfi_: p; ] AS 

1 e 2 

- fj - i ax - a 6 e 2 i /e x 



\ 


■SAP 

d X 


IaL « [Ax - Afej/fj 

d x 


( 6 - 2 . 8 ) 


For the astronomical observations, it is necessary to re-define, 
slightly, certain vectors in the table of vector definitions on Page 64-C 
Namely, the vector is computed as 


and is the position vector of the target body relative to trucking 
Station 1. Station 1 will now be taken to be the astronomical observ- 
atory from which the right ascension ( ^ ) and declination ( £ ) 
measurements are mad' In order to compute the apparent ol and & , 

P l must rotated from its ecliptic representation into the geo- 
equatoribi : ?ccoHing to the transf ormation 


* t 


From the components of ^ , the right ascension ^ target body is 


computed ai 


- tan" 1 ( Ply / ^lx) 


and the declination as 


iin " 1 ( flz /t\ 


where ^1 * I Si l . 
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The observation partials for the astronomical observations can be 
written as 


H >«,S) 1 

-f.l (<_._*> 1 



[ *(JL, . *p>. 

(2x6) [ 3 <-; , ^ ) 


t 3(1 . V 


where the :ond ter*" on the right hand side of this equation is identified 
as 






1 


Elements in the first matrix on the right hand side above are given by 


3* = - sinOt / ( cos £ ) 

3 x P ' 


3 * = 

3y; 


+> cos o L / ( cos S ) 



i 

» j - 
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9«l - 9cC . 3*. « - 0 

d xj» a*J T5J” 3z£ 


3 S = - cos*/, sin S / (*. 

TxT 1 


11 = - sinet sin $ / P 1 

ay/ 


a s = cos s / e, 

a 2 ; 


a s = a s = *v s = o 

a*; *y air 


These partials are handled differently than the observacicn partials 
for other deep space data types. Since the astronomical observations are 
made on the target body whose state must be augmented to the usual S/C 
state vector as solve-for or dynamic consider parameters, the partials 
correspond to partitions H g or H^, respectively, in the augmented 
observation matrix. Furthermore, only when there are onboard optical 
observations or significant dynamical interactions between the S/C and 
the body is the S/C knowledge covariance affected by astronomical 
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observations. One final note about astronomical observations: G0DSEP's 

mathematical model has been designed to handle both measurement noises and 

biases separately for each observable ( oL and/or $ ). 

The following definitions are used in the azimuth and elevation 

angle partials. 


$ 

e 


A 

e 


§p 


S/C azimuth, measured positive from north toward east 
(See Figure 6-2) 

S/C elevation 

S/C range vector from station 

Unit vector in direction 

Projection of £ onto plane normal to x 


= Geocentric equatorial S/C position 


= Heliocentric ecliptic S/C position 


x 

—s 


Geocentric equatorial station position 


Unit vector in x direction 
— s 


A 

w 


A 

u 


Unit vector orthogonal to x g and pole (local east 


from station) 

* Unit vector orthogonal to x^ and w (local north from 
station) 

l = Transformation from equatorial to ecliptic coordinates* 

For simplicity, all azimuth and elevation partials are derived in 
geocentric equatorial cartesian coordinates and then transformed to 


ecliptic. 
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Again referring to Figure 6-2, the projection of £ onto $ will have 
magnitude cos | sinol , or 

sin< *= sec Q [ a 1 T p 


where 


a* 

tan ^ tan ^ 
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+ 

v *' 
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(6-2.14) 


For use in Equations 6-2.11 and 6-^. 13 above 
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(6-2.15) 


where x^ and x ? are components of x^* Finally, the partials 
s s 

computed in equatorial coordinates must be transformed into ecliptic 


2S+JLL - P ) . 
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( 6 - 12 . 16 ) 


For vehicle based optical measurements, we use the following 
definitions. 

x • spacecraft helio 9 entric cartesian position 




plane '/tar get body heliocentric cartesian position 


• x 1 planet range vector 




d “ vector of star direction cosines 
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S 

T 

$ 


vector of target planet orbital elements (Keplerian) 
target planet radius 

star-planet angle 

apparent planet diameter measurement - angle subtended 
by planet disc at the spacecraft 

re ro vector, 3x1. 


Star-planet angle pertials: 


cos 


a «r £ 

d T e 


( 
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\t Apt ' 5 


- cos 


v e ) 


(6-12.17) 
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and if the ephemeris body elements are Keplerian rather than Cartesian 
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where the partials of cartesian to Keplerian elements are computed 
numerically. 

Apparent planet diameter partials: 

.i. i * V P 
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For any data type which has a bias, 


y ■ H £ X + b 



1.0 


( 6 - 12 . 20 ) 


( 6 - 12 . 21 ) 


( 6 - 12 . 22 ) 


(6-12*23) 


6.4 FILTER 

After the knowledge covariance h*s been propagated to a measurement 
time, and the observation matrices are computed, the 0D filter can per- 
form its function of estimating the set of solve-for parameters (non- 
deterministically) and updating the knowledge covariance accordingly. 

There are two types of filters available, Kalman-Schxnidt (K-S) and 
weighted least squares (WLS), pluc capability for a third filter, to be 
established by the user. K-S is the most commonly used filter because it 
treats consider parameters in a realistic fashion. 

The filter updating process requires computation of several matrices. 
i First, the propagated estimation error (knowledge) covariance F at the 

measurement event time is 
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x 

: T 

xs 


xs 


xw 


XV 


P c 

V vw 

C T P 
vw w 


where P , •••, P are the covariances of the S/C state, ignore para- 

x w 

meters, respectively, and C ,•••,€ are the cross-covariances between 


xs 


vw 


appropriate augmented parameters. The observation matrix H defined in the 
previous section is 

‘ H 


H 


H 

s 

H 

u 

H 


l w 


The measurement residual matrix J is defined "as 

J - HPH T + R 

where R is a diagonal matrix containing variances of the measurement white 
noise. For example, a simultaneous 2-way/3-way range measurement would 
look like 


R - 


** 2 


*2R + S’ 


J 3R 


2 x* 2 

where 0^ Is the 2-way range noise variance and 0^ is the additional 

3-vay range noise variance due to timing synchronization. For a single 

star-planet angle measurement, R would be a scaler 

R - 0* + 0j?/r 2 
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2 2 

where end (5^ are the optical resolution and planet or body 

center finding noise variances, respectively, and r is the spacecraft 
range to the planet. 

The updated covariance (P + ) after the measurement has been 
processed (and in theory after the sta^e estimate has been updated) 
is in general 


P + - (I - KH) P (I - KH) T + KRK T 


(6-3) 


where K is the iilter gain. 


KALKAN-SrHMIDT 

The filter gain for K-S is straightiorward 
K » PH T j' 1 


Since only estimated parameters can be updated by the OD process, 
the entries of K corresponding to consider and ignore parameters 
•re zeroed out, that is, K « [ K x K g 0 0 0 ] T . The updated 
covariance is then formed by Equation 6-3. 



WEIGHTED LEAST SQUARES 

The sequential, or recursive weighted least squares (WLS) algorithm 
implemented in GODSEP is equivalent to a batch WLS filter if there is no 
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process noise. Since process noise is a significant part of low thrust 
analysis, the WLS filter muot be used recursively, because it has no 
batch equivalent. The sequential WLS consider filter acknowledges con- 
sider parameters only in the covariance update (Eqn. 6-3) and not in the 
gain matrix calculation. Therefore a set of reference covariances 
for the state and solve-for parameters must be maintained at all times. 

This sec also represents the filter analysis as it would be in non-consider 
form. 


Thus, the WLS filter computation requires three operations: (1) pro- 

pagation of the reference covariance to the measurement event, (2) com- 
putation of the filter gain and (3) updating both the reference and 
knowledge covariances. 

(1) The reference covariance (P) consists of 


A 

P 


xs 


xs s 

and is initialized at the a-priori values. Thereafter, it is 
propagated from one measurement to the next by 


where 4 is the augmented transition matrix corresponding to the 
x and s_ parameters. P is computed in parallel with the actual 
knowledge covariance P. 

(2) Given $ at the measurement event, the WLS filter gain is 


A 

where H 


ft] 


A A T A « »T -1 

P H (H P H + R) 
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(3) The reference covariance is updated, after measurement pro- 
cessing, by 

P - (I-KH) P 

and the knowledge covariance P is updated by Eqn. 6-3. 

Thus, at measurement events, the 0D filter updates the knowledge 
covariance to simulate taking a tracking measurement, processing the measure- 
ment in an orbit determination algorithm, estimating desired parameters and 
reducing (or updating) the knowledge uncertainty to reflect this new infor- 
mation about the trajectory. 

6.5 GEKEt&LIZED COVARIANCE 

The function of any filtering algorithm is to process available measure- 
ment information and produce a best estimate of the spacecraft state and 
any parameters that are being solved-for. Best is usually defined in a sta- 
tistical sense, such as the minimum variance processes used in differing forms 
in the weighted least squares and Kalman-Schmidt filters. But in practice, 
filter performance is dependent on how well the assumptions uned in the 
filter definition approximate real-world processes, because all error sources 
cannot be modeled, nor can those that are modeled ever be modeled exactly. 
Therefore, each filter must be evaluated not only on its ability to produce 
small error covariances in the resulting estimated state, but also be as 
insensitive as possible to errors in its model assumptions. 

Generalized covariance error analysis is a useful tool for studying 
filter sensitivity. For generalized covariance studies, two sets of know- 
ledge errors are carried during the orbit-determination process. Assumed 
knowledge uncertainties are those generated by the filtering algorithm 



according to the mathematical model and all the assumed errors input to 

* 

It. True knowledge uncertainties represent the effect the filtering 
algorithm has on actual state estimation when the real-world error sources 
are not the same as those assumed by the filter. Evaluating filter sensi- 
tivity to a model assumption involves comparing the resultant effect on 
assumed and true uncertainties of a modeling mismatch between the filter 
and real-world uncertainties. This modeling mismatch is accomplished in 
one of two ways; true a priori uncertainties may be set at levels other 
than assumed levels or the true state may be augmented by a vector of ignore 
parametcrs--parameters whose uncertainties arc recognized by the true 
covariance analysis, but which are completely ignored by the assumed filter 
ana lysis. 

The filter that is least sensitive to a model mismatch is determined 
on the basis of two criteria. First, which filter yields the smallest true 
estimation errors. Second, for which filter are the true errors most closely 
approximated by the errors predicted by the 'filter covariance analysis. Thus, 
given a mission with a specific set of model mismatches, if two different 
filters produce equivalent true errors, then the superior filter is the 
one whose assumed errors are closest to the true ones. Similarly, if the 
resultant assumed errors of the two filters arc equivalent, the superior 
filter is the one with the least true estimation uncertainty. Generally, 
qualitative judgments are required because several sets of mismatches must 
he studied, and the relative performance of the filters may vary. 

In error analysis, generalized covariance is a filter sensitivity 
study tool that is normally available only in a simulation program. It is 



accomplished in G0DSEP with a minor increase in core and comoutatlonal 
time compared to a full simulation and has the additional advantage of 
generating ensemble true state statistics rather than a single sample as 
in a simulation. The only disadvantage of generalized covariance i» that 
it uses the same linearized dynamic and observation models as the assumed 
filter analysis, and can therefore not study problems that arise from 
nonlinearities. 

The actual operation of generalized covariance in G0DSEP requires tha 
a standard error analysis be run first. The filter gains, associated with 
the assumed knowledge uncertainties, are stored on disc or tape. Now the 
error enalysis with all the same measurement events is repeated. Only thl 
time, a-priori uncertainty levels and measurement noise are modified, and 
ignore parameters are added, to the extent of desired mismodeling. At 
each measurement event, Eqn. 6-3 is applied to what is now the true know- 
ledge covariance using the appropriate stored filter gain. The true 
covariance analysis thus proceeds in analogous fashion to the previous 
assumed covariance analysis. Obviously, many mismodeling conditions can 
oe studied with the same filter by repeating the generalized covariance 
analysis. 

6.6 GUIDANCE 

Although the knowledge covariance is modified by measurement events, 
the control covariance, which represents the ensemble of acrual deviations 
from the desired or reference trajectory, will grew without bound. The 
only process which will reduce the control covariance is a guidance event 
that is, the design and execution of a trajectory correction, either 
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impulsive AV or low thrust . 

Low thrust guidance represents an update of the nominal thrust con* 
trols (magnitude, direction and cut-off time). In terms of system cost 
and efficiency, it is better to use the existing low thrust propulsion 
system for guidance than to add auxiliary means, for example high thrust 
chemical engines to produce impulsive &v. Of course, certain problems 
inherent in low thrust propulsion, in particular terminal controllability, 
may force the addition and use of an auxiliary chemical propulsion system. 

In mathematical terms, given a trajectory state deviation, gx^ ® 

S 21 (t Q )> where t* is the guidance epoch, we wish to null out the effects 

of S by making a bias type correction $u to the nominal thrust controls. 

To be efficient, the correction is applied over some finite interval 

[ t , t I such that the target error £t , caused by Sx , at some final 
o c J *— ‘O 

time (tj) is removed. For linear analysis, we seek the guidance matrix T , 
such that 

sb - r ^ 


The linear ensemble of thrust control corrections is then 


U - E [gu Su T ] 

u « r sx/] r T 


(6-4) 


In G0DSEP, the trajectory error ensemble E £ S $x q Jis the control 
covariance P c (t Q ). Using P^ represents a pessimistic sizing of the 
thrust corrections because only the known trajectory error (not to be 
confused with the knowledge error) can be removed. The known error 
generally corresponds to the control error as long as P c ”(t o )5^ P^(t Q ) 
where P^ is the knowledge covariance. 
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To compute the guidance matrix 
matrices 


P we first compute the sensitivity 


S(t f ,t # ) 


9 1 < tf ) 


9 


U 


a i (t f > 

» * <t G > 


or. 


S 


V 


f jli V J 
L 3x(t f ) J 

r 1 

l J 


♦ (t fV 9 

♦ <t f ,t c ) t <t c ,t 0 > 


where the first matrix in S and V is formed by numerical differencing 
and the second two matrices (<^ and 9) sre obtained from transition matrices 
generated by the trajectory propagation routine (Section 6.2). If variable 
time of arrival is desired, the control array $ u is augmented with the 
arrival time and the S matrix is augmented by x (t^), relative to 
the target. 

The guidance matrix can now be defined by 
P - -W^ 2 S T | sw u 2 S T ] V 

P - -[s T w 2 s] s T 


IfN u >N T (6-5.1) 

If N t > N u (6-5.2) 


where and are the number of parameters in the control and target set, 
respectively, and and are diagonal weighting matrices for the control 
ar.d target parameters, respectively. The first f^rm of T, 6-5.1, reflects 
a tnir mum quadratic control correction and the second, 6-5.2, corresponds 
to minimum quadratic target error. 
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Generally, there are more control parameters than targets with the 
exception of two cases: (1) terminal approach to the target where con- 

trollability drops off rapidly; and (2) the application of control 
constraints which effectively reduces the set of available controls. 

If constraints on the control corrections are imposed, then P mast 
be suitably modified. First, the unconstrained P is computed along with 
the ensemble unconstrained control corrections, 


U = p P c " P 


( 6 - 4 ) again 


2 

Each liagonal component of U, , is compered against its constraint 

value ( . If 0* u is greater than £ u max’ then the a PP ro P riate 

row of T is scaled by Su„. v /<Tu. The total control set (and guidance 
matrix) is then separated into two subsets: unconstrained controls, & u j 

and P^, and constrained controls, £ an< * ^ 2 * 

r, 


= c 

% Ho 


max r 

*u ' 2 J 


The new control corrections are computed with (6-4), (actually only the 

remaining unconstrained controls are computed ) the test for constraints 

are made, T id modified again, and the entire process is repeated until 

all constraints are met, or there are no more controls left. The guidance 

corrections are executed (figuratively) at time t^, that is, are uplinked 

to the spacecraft, but apply over the entire guidance interval [ t Q , t c ] . 

Execution of the guidance updates causes the control covariance to 

diminish from t to t whereupon it begins to grow again. Guidance accuracy 
o c 

is measured by how much the control covariance can be reduced at t £ which 
depends upon how well the thrust corrections were designed. The limiting 
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accuracy of maneuver design is the knowledge error or covariance at t Q . 

Thus, the post maneuver control covariance at t £ is the propagated 

knowledge covariance (as per Eqn. 6-1) from t to t , 

o c 

P + ( k ) - f(t ,t ) P, (t ) f (t ,t ) + Q(t ,t ) (6-6) 

c v c r c * o' k o » v c* o c* o 

In G0DSEP, to denote guidance execution, P*(t ) is set equal to 
P^(t o >* This is equivalent in effect at t c to applying Eqn. 6-6. 

However, it means the value of P+ in the guidance interval is not valid. 

This is a relatively minor problem compared to the reduced burden on 
computational storage and logic. 

One exception to setting P+ = P^ occurs when there are more 
targets than available controls, which often happens when control 
constraints have been activated. In this case, there will be some non- 
zero target error that was not removed by the guidance corrections. This 
implies that not all of P^ was removed. Hence, the post maneuver control 
covariance must include the residual state error. 

p c = p k + { vT t wT 3 l v + s r] } p c ' { v T [ w T j I v + s dj 7 

As long as no more guidance events are executed between [t Q ,t c J , 
updating the control covariance at t is theorettically no different 
then using Eqn. 6-6 at t c « However, if another guidance event is scheduled 
in the burn interval, say at t^, then a somewhat different logic is applied 
to size the correction. It is assumed that the first guidance 'vent between 
t Q and t is a primary maneuver. Subsequent guidance events in this inter- 
val are considered to be vernier maneuvers, representing refinements of the 


thrust corrections computed in the primary maneuvers. 



For vernier maneuvers, 9 and ^ are redefined for the vernier burn 
interval and P is computtl as in Eqn, 6-5, The guidance updates are 
computed using Eqn. 6-4 with the knowledge gained since the primary 
(or previous vernier), that is, 

E [*Sl *41- <<'l> - W 

Recall from the previous discussion that P+(t^) * s usually the propagated 
knowledge from the previous guidance event. 

A measure of guidance effectiveness is the estimated target error 
before and after the maneuvers. 

before guidance correction, E t St t ] * vp' (t Q ) v T 

St St J = v p (t Q ) v . 

This simple measure assumes, of course, that no further dynamic error 
* will occur from t to the target time t . 
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An Important part of the guidance ar.d navigation process is 
the time interval from the last navigation measurement used for 
guidance design to the actual time of guidance implementation 
(t o ). The time interval (or delay) is necessary for ground pro- 
cessing of all previous measurements, estimation of the S/C state, 
designing the thrust updates to correct the trajectory, and 
execution of the updates. Typical intervals are 3 to 13 hours, 
and are usually critical only in the terminal mission phase where 
trajectory controllability (with respect to thrust controls) 
diminishes rapidly. This time delay is user specified for each 
guidance event. 


Impulsive AV guidance is very similar to low thrust except 

for a zero burn interval (t « t ). The deltas velocity is treated 

o c 

as if it were a control correction £ u, that is, 

a» - r Si. 


To compute P , the sensitivity matrix S is first partitioned 

into position and velocity submatrices, 

« 

S -[*:»! 


J T(t f ) 


ai(t £ ) 


where A 


and B « 
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If the target vector T has only two components, e.g., B*T and B*R, 
then 1 | is minimized and 

V 

T - [ -b t <bb t ) a j -b t (bb t ) b ] 

If T has three components, e.g., the position vector at t^, then 

r - 

The computation of ensemble velocity correction, U in Equation 6-4, 

follows directly. As in low thrust guidance, the pre-maneuver control 

covariance P (t ) is used to size U. 
c o 

i - e[a* av t ] - p ?; <t o ) r T 


Execution errors related to low thrust control updates are neg- 
lected because they are second order effects compared to thrust error 
associated with the nominal thrust profile. However, AY execution 
errors are taken into account because impulsive maneuvers often occur 
during ballistic or coasting portions of the mission and can represent 
a significant contribution to trajectory error. In order to compute 
execution errors, the most probable AY %s first determined by 
the Hoffman- Young approximation (Reference 8). Let A A 2 » Aj 
be the eigenvalues of the AV covariance, U, and be the largest 
eigenvector of U, define 

A - >j+A 2 + A 3 

- fs: f l+ M™ll 

'f T I J 
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I 



e 



I 


82 


then the probable AV Is e[ay! 


0 & 


AV X 
AV 2 
1 ^3 


Now the 3 x 3 AV execution error covariance Q is composed of 


* - AV, 2 (T 2 + AV 2 J p 2 o; 2 + ^V, 2 av 3 2 <r, 2 

11 1 m cT 1 

“xy 


(6-7) 
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where P * 

4 xy 

<T . 2 






AV resolution and proportionality variances 



ecliptic (X-Y) pointing variance 


out of ecliptic (z) pointing variance. 


As in low thrust guidance, the post-maneuver control covariance 
(t Q ) is set equal to the knowledge covariance P^ (t Q ) corrupted by 
the AY execution errors, 


K <«„> 


p. (t ) + 
k o 



Pre and post maneuver target uncertainties are computed in equivalent 
fashion to low thrust guidance. 



7.0 SIMSEP Analysis 


The trajectory simulation mode SIMSEP has been designed to 

* ^ 

provide deterministic analysis of ballistic and low thrust missions. 
Computationally, SIMSEP imitates "real" trajectories in the presence 
o£ a wide variety of environmental and system uncertainties. A pri- 
mary objective is to deduce expected or probabalistic behavior of 
the real mission by studying a relatively small subset of simulated 
missions. 

The purpose of this section is to discuss the key analytic 
concepts in SIMSEP. This will be done in two parts: 1) by discuss- 

ing the principal algorithms, and 2) by outlining the basic compu- 
tational structure. Although many algorithms used in SIMSEP are 
similar in function to algorithms used in T0PSEP and C0DSEP, their 
specific'applications here warrant an extended discussion of their 
underlying theories. 

7.1 Program Scope and Methods 

Before proceeding with a step-by-step description of the 
algorithms and computational structure, it is worthwhile to compare 
the essential similarities and differences between G0DSEP end SIMSEP. 
Unlike the error analysis mode which works exclusively with error 
ensembles and a reference trajectory, the simulation mode actually 
formulates many discrete examples of the "real world" or "actual" 
trajectory. Each of these is propagated, In a deterministic sense, 
by the sat^e trajectory integrator, integrating the same equations of 
motion. However, many variables and parameters appearing in these 
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equations are subjected to random alterations, corresponding to 
discrete uncertainties. Hence, each deterministic simulation f a 
mission is different according to the effects of the sampled errors. 

Operationally, SIMSEP does not sample each error In succession 
and propagate trajectories with just one error source active at a 
time. Rather, all are initialized by differing amounts for each 
actual trajectory. Thus the averaged effect of all error sources 
acting in concert can be estimated by repeating the mission simula- 
tion process a sufficient number of times. This is the essence of 
the Monte Carlo method and is the basis of the simulation approach 
to determining how trajectory nonlinearities and uncertainties can 
affect the G&N process. 

Perhaps the best example to illustrate the fundamental differ- 
ences between the methods used in G0DSEP and SIMSEP is the problem 
of propagating, or mapping, an error covariance from one point to 
another along the reference trajectory. It-will e recalled that 
the principal method for propagating a covariance in G0DSEP is by 
the state transition matrix mapping, namely, 

P k+1 “ fk+1, k P k tlefl, k 

where represents the state transition matrix and P^ and 

P^j are covariances at t^ and tj^, respectively. When there is 
dynamic process noise, an effective process noise matrix, ^ 

la also added; (See Section 6.2). The state transition matrix is 
generally computed simultaneously with the trajectory by integrating 




92 


the variational equations. However, the variational equations are 

based on linearized expansions of the differential equations of 

motion and neglect all second and higher order terms. Hence, a 

% 

state transition matrix mapping of covariances must theoretically 
be limited to mapping covariances within the envelope of linearity 
surrounding the reference trajectory. Rarely is this assumption 
true at all points along an Interplanetary trajectory, especially a 
low thrust trajectory. Covariances propagated by this means are 
subject to error whenever a region of significant trajectory non- 
llnearltles is encountered. 

On the other hand, the method for propagating a control 
covariance in SIMSEP is not plagued by these effects, although it 
has its own peculiar shortcomings. The SIMSEP approach to this 
problem relies on the Monte Carlo method where a multitude of 
sample trajectories are propagated between the two time points in 
question. The trajectory state vector data at -t^+i are processed 
and accumulated in such a way that the covariance can be reconstructed 
by standard ststistical calculations. Hence, SIMSEP maps a covariance 
as a statistical ensemble calculated from many data points and not as 
a simple mathematical entity like G0DSEP. 

Therein lies the primary drawback to the Monte Carlo method and 
to the use of SIMSEP for general G&N analysis. Although the Monte 
Carlo method will, in theory, converge to the exact covariance, the 
rate of convergence tends to be extremely slow. For accurate 
statistics, inordinately large amounts of computer time are often 
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necessary to perform the many trajectory propagations. 

Although both G0DSEP and SIMSEP are nreflight mission analysis 
programs used to identify the general G&H subsystem characteristics, 
they are usually used at diff'-ing phases in the development of an 
overall systems analysis. For the purposes of a preliminary systems 
design, a G0DSEP analysis is the most cost effective means of evaluat- 
ing the basic G&N subsystem requirements. As such, SIMSEP is generally 
relegated to verifying the linear analysis results, and only in the 
advent of serious nonlinearities is the simulation mode called upon 
for more extensive studies. 

7.2 Definitions and Concepts 

The first important concept in SIMSEP and common to all MAPSEP 
modes is the reference trajectory, denoted by » 
reference trajectory is computed under some set of "reference 
integrating conditions" to satisfy desired targets at mission's end. 
Moreover, represents a deterministic solution to the equations of 
motion for the assumed dynamic and systems models. For SIMSEP, the 
Initial state and reference integrating conditions, i.e., ephemerls 
parameters, thruster characteristics, etc., are read as input since 
it is assumed that they have already been computed as output from a 
T0PSEP analysis. 

A second quantity important in SIMSEP and common to G0DSEP is 
the control erroi covariance, P^. Generally, an a priori control 
covariance is defined at Injection (or at the starting point of the 
mission being studied). This siatrlx mathematically describes the 
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distribution of real state errors relative to the initial reference 
trajectory state. In SIMSEP, it is implicitly assumed that the 
probability distribution of these errors is Gaussian with zero mean. 
Once an a priori control has been given, it is randomly sampled to 
form an error vector, S X^» which corresponds to a deviation of the 
actual trajectory state relative to the reference. At the same 
time, error sources associated with the host of other dynamical 
and systems uncertainties are also sampled to create the so-called 
"real world integrating conditions". For the actual trajectory 
state vector, this procedure may be written as 

h ' * r + *h 

where $2^ is a deviation obtained by sampling Utilizing 

these integrating conditions, the actual trajectory state, 3^, is 
propagated from point t^ point as a discrete example of an actual 
trajectory. 

The third-ycritical variable used in both GODSEP and SIMSEP is 
the knowledge error covariance, This matrix is propagated in 

the error analysis from measurement to measurement where it is 
systematically updated according to the filtering algorithm. In 
SIMSEP, instantaneous evaluations of P^ are input at each guidance 
event and are left unchanged throughout a given run since there is 
no explicit orbit determination process modeled. However, the 
knowledge covariance, like the control covariance, is sampled to 
formulate an error vector. This error vector determines the error 
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in the estimated state relative to the actual trajectory state, and 
is used to compute an estimated state vector by > 

h m h + z 

where e is the sampled error vector from P^. If other parameters 
are estimated during the orbit determination calculations, they are 
included as augmentation parameters. These too are sampled in order 
to formulate a set of "estimated world integrating conditions". 

With each of these key quantities having' been defined, it is 
worthwhile to mention why each is important in a simulation run. 

The reference trajectory, for example, serves to define the mean 
for all actual trajectories, as well as defining the reference target 
conditions used during guidance. The actual trajectory is, of course, 
the mathematical representation of the real motion and is carried 
from event to event until the final target is reached. On the 
other hand, the estimated trajectory is used - exclusively for re- 
targeting the actual trajectory back to the desired targets and is 
computed only during guidance. 

7.3 Guidance 

One of the principal purposes of SIMSEP is the detailed exami- 
nation of nonlinear trajectory effects, especially as they bear 
upon the guidance problem. In this section, the fundamental concepts 
underlying linear and nonlinear guidance will be presented, and the 
implementation of these concepts into algorithms will be discussed. 
Beforehand, a careful distinction between targeting and guidance 
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must be drawn since MAPSEP has algorithms for both and since many 
of the basic steps and operations are similar. Whereas a targeting 
problem is solved by formulating an entire control strategy for a 
complete mission, a guidance problem assumes that a solution to the 
targeting problem has already been found. Furthermore, the current 
trajectory which is to be corrected is assumed to be in a "close 
neighborhood" to the original reference solution. Hence, the con- 
trol changes computed by a guidance law are expected to be small 
refinements to the original controls, even in the presence of non- 
linearities. 

7.3.1 Linear Guidance 

The linear guidance option in SIMSEP is analogous to the guid- 
ance used in G0DSEP except that it applies to a discrete trajectory 
error as opposed to an ensemble of errors. For both modes, the 
linear guidance matrix is the same. To compute a guidance matrix, 
a sensitivity matrix is evaluated between the point of the guidance 
event and the target about the reference trajectory. This matrix 
of linear partials relates control changes to target deviations and 
is used to map estimated trajectory errors, into control 

updates, according to the guidance laws: 

1) A V * P 6Xg , for impulsive corrections, and 

2) Au ■ $Xg , for low thrust. 


In spite of its overall simplicity and ease of implementation. 
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the advantages of linear guidance are off-set somewhat, by its draw- 

i 

backs. First, the trajectory error, S X^ . must lie within the 
envelope of linearity for this to be a valid method. Whenever this 
is violated, the resulting guidance correction can be invalidated. 
Furthermore, a linear guidance correction is executed without itera- 
tions. In fact, with this guidance there is no direct evaluation of 
a control correction's effectiveness in reducing target error. Only 
if the updated trajectory is propagated to the target c.n the result- 
ing target error be determined, and even then there is no recourse 
for making further corrections if the original correction is 
ineffective. 

7.3.2 Linear Impulsive Guidance 

The essence of impulsive guidance is founded on the mapping 
relations which propagate arbitrary linear deviations relative to 
some known trajectory into new deviations at some later time. 

Clearly, this is a property of the state transition matrix which 
maps a six component state deviation, £x^ . evaluated at t^ into a 
new deviation, S X,^ [ . at t^+j, by the equation, 

•^k+1 “ $k+l,k -^k 

If is the target time and t^ is the time of the guidance event, 

then k can also map state vector changes, like an impulsive 

velocity correction, into state changes at the target. 

However, in most analysis the actual target conditions are 
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specified in terms of target variables such as B-plane parameters, 
Keplerian elements, etc., instead of X, Y,...., Z state coordinates. 
Fortunately, the target variables are functions of the final trajec- 
tory state and it is possible to generate a differential transforma- 
tion of the form 


*r - ^ ix kll (7-2) 

which transforms differential coordinate changes into target variable 
variations. In the above equation, represents a matrix of linear 
partials of the form 


_ j) Of T 2 T n> 

V * (x, y, z, *, Y, 2)^ 

where there are n-target variables, T^, T£ and six state 

components. By substituting Eq. 7-1 into 7-2, a relation for mapping 
state changes at t^ into target changes at t^ + j is obtained, that is, 

“ *1 Jk+l.k -^k “ 


Performing 

fit+i.k 


the indicated matrix multiplication and replacing 
with N, the equation becomes 



(7-4) 


where N has dimensions (n x 6). 

With this background, the impulsive guidance problem can be 
stated most simply as a determination of a velocity change, 4V. 
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which) when added t< Me- nulls the target error, &T . Again 
writing Eq. 7-4 but in a partitioned format, 

»- 

St - N(l) + N(2) Sy k , 

it is recognized that % T can be made zero by adding some approprl* 
ate AV to *••• 

N(l) $ r fc + N(2) (tv k + £0 - 0. 


For the case of three-variable impulsive guidance, i.e., three 
unique targets, the solution for AV is given as 


AV - - 1U2)' 1 N(l) $r, 


provided N(2) is nonsingular. Note that this can be re-written as 


4 H - ^- N < 2 >" 1 NCD i - 1 3 


4/ - r u. 


<7-5a) 


(7-5b) 


the desired guidance law. 

For the case where there are two target variables Instead of 
three, the problem has more controls (3-velocity components) than 
end conditions and a generalized inverse which minimizes the magni- 
tude of the velocity correction is used according to 



* . « 11 - 



i. I) 
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- N(2) T [n( 2) N(2) T *J 1 N(1) -^k 


- m'Y i£* 


where N(l) and N(2) are non-square matrices with dimensions (nx3). 
Again this relation can be re-written as 



(7-6a) 


(7-6b) 


Algorithms based on Eqs. 7-5a and 7-6a are the basis of the 

linear Impulsive guidance contained in subroutine LGU1D. The guid- 

•n 

ance matrix, \ , for either the two or three variable cases, are 

computed as outlined above and the state vector deviation, Sjl. 

Is calculated as the error in the estimated trajectory state relative 
to the reference, namely, 

Aik * * h - h’ 

evaluated at the guidance event. 

7.3.3 Low Thrust Linear Guidance 

The low thrust linear guidance law has the same format as the 
impulsive law except control changes are made to the vector of low 
thrust control variables, *11 . Another difference is that the low 

N 
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thrust acceleration acts slowly to bring about state changes; hence, 
the low thrust linear guidance ruatrix operates in an integrated 
fashion over a fixed tralectory segment to redirect the motion. 
Otherwise, low thrust guidance is simply an extension of the basic 
methods discussed above. 

The most complex part of computing low thrust corrections is 
the determination of the guidance matrix, T* 1 • This matrix depends 

not only on the trajectory dynamics between the maneuver point and 
the target, but it also depends on the trajectory response to con- 
trol changes, i.e., controllability. As before, the first step is 
to integrate the reference trajectory from the guidance point to 
the target, evaluating the augmented state transition matrix. I* 
SIMSEP, the transition matrix is computed by integrating the vari- 
ational equations as was discussed in Section 4-5. By selectively 
partitioning the transition matrix, the requisite sensitivity matrices 
relating state and control variable deviations to future state devia- 
tions are obtained. 

In terms of partitions in the augmented state transition matrix, 
state deviations at the target time, are given by 

-^k+l " $k+l,k -i^k + ® ^ ^ 7 " 7 * 

where | k is a state transition matrix as defined in Eq. 7-1 
and 0 is a matrix which maps control variable deviations into 
state changes at t^ + ^. S u in Eq. 7-7 corresponds to a set of 
thrust control biases. © can also be written as 



'KlWHNpjw#. wiiu w»>Wh<w> 
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i > 


0 - ^ (x » Y> z> *> 1 > *> k +i 


(7-7.5) 


and Is seen to be (6xm) where m is the number of controls. Follow- 
ing the same line of reasoning as was given in Section 7.3.2, it is 
recognized that partials of target variable variations with respect 
to control variable changes are needed. Hence, Eq. 7-7 is multiplied 
by the transformation matrix (Se 7-3) to obtain, 

St = ^ f k+l k +-*1® ii. 0-8) 

Therefore, the guidance problem is reduced to finding a /\u w hioh 
when added to &U will make St = 0. For convenience, it is assumed 
that & u is either zero or that it can be solved-for during the orbit 
determination; thus permitting Eq. 7-8 to be re-written as 




© Aik + /v » \ 


k+1 , k 


SX. = 0. (7-9) 


For the problem where the number of controls (m) equals the number 
of targets (n) and the matrix ^ has an inverse, the solution for 
A u in Eq. 7-9 is 

-1 


^ m ® !k+i,k 


(7-10) 


( 1 


where £-*e is the deviation of the estimated state vector relative 
to the reference at the guidance event. From Eq, 7-10 it is easy 
to see that the desired guidance matrix for this particular case 
is given as 



i- . . 


1 
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V 

and Is a (6x6) matrix. 

The problem is complicated somewhat when the number o£ targets 
la less than the number of controls. In this case, a generalized, 
or pseudo-, inverse matrix operation is used, and the transformation 
matrix does not drop out. Nevertheless, a solution is obtained 
by determining the Au. that makes St * 0. Letting A "’’l® and 
B “ ^ k in Eq. 7-9, a particular solution (out of the 

infinity of solutions) is given as 

Ah - - a t £aa t ] b $Xg. (7-11) 

This particular choice of AH also minimizes the magnitude of the 
control change (See Section S.3). Therefore, the desired guidance 
law can be written as 

a* - r 

% • 

where T * • t AA* 3 B. 

As before, the computational steps described here are imple- 
mented in LGUID. 

7.3.4 Nonlinear Guidance 

Nonlinear guidance parallels in many respects a targeting 
problem where an iterative, linear algorithm is used to determine 
control changes. The primary difference is that the trajectory is 
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assumed to be reasonably close to the reference. By reasonable, it 

* 

Is meant to imply that the algorithm should be able to redirect the 

/ 

s/c motion to the designated target by making a few Iterations (three 
or four). In practice, a real trajectory can deviate widely from 
the reference and thereby require as many as eight to ten iterations 
before the guidance algorithm is able to compensate. For situations 
where convergence is not achieved after many iterations, the guidance 
Is said to be divergent. The real criterion for qualifying a maneuver 
as divergent is somewhat subjective and established by the analyst. 

In some cases, an extremely slow rate of convergence on the part of 
the linear correction scheme can be attributed to non-adaptive itera- 
tion logic. 

'Mathematically, divergence implies that the real world integrat- 
ing conditions acted in such a way that the guidance algorithm was 
unable to rectify the motion. Physically, divergence suggests that 
something in the dynamics or s/c system has been mismodeled, or under- 
designed, snd that it has interactions with the other systems to cause 
vide, usually nonlinear, deviations from the reference mission. From 
the G&N point of view, it is these missions that are often of the 
greatest interest. They identify potential problems either in the 
baseline configuration or with the navigation and operational proce- 
dures. In many Instances, divergence in the guidance can be traced 
to some well-understood phenomena, e.g., controllability, trajectory 
non-linearities, suboptlmal schedule of guidance events, etc., but 
a clear identification and resolution of these problems in terms of 
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changes to the system design and/or operational procedures may not 
be so straightforward. 

The basic computational steps taken in a single iteration are 
as follows: First, an estimate of the actual state vector and the 

corresponding estimated integrating conditions are obtained from 
the simulated orbit determination logic. The estimated state is 
integrated to generate the estimated trajectory between the guid- 
ance point and the target. At the targeted stopping conditions, 
an estimated target error is computed by 

A I = Ir - Tg 

where T_ and are the target variables on the estimated and refer- 
ence trajectories, respectively. Next, a sensitivity matrix, S, of 
target variations with respect to control variations is computed 
about the estimated trajectory according to the matrix relation, 

S - T e 

where *} is the state to target transformation defined in Eq. 7-3 
and where 6 is that partition in the augmented state transition 
matrix which maps control changes into state variations (Eq. 7-X5). 
The matrix $ has the format, 


s 


3 (T, T 

Li- ii 

3 u 2» 


T ) 


t i±jj- 


V 

m 


and has dimensions (nxm) , where n is the number of targets and m 
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4 

the number of controls. With S it Is possible to relate control 

♦ 

changes at the maneuver to target variable changes according to 

t 1 - S A u . (7-12) 

If the matrix 5 is square, i.e., the number of controls and 
targets are equal, then the solution to (7-12) r ;>'res only a 
single matrix inversion, namely, 

& JL " S' 1 AT. <7-23) 

However, in most practical cases, S a nonsquare matrix (m>n) 

and a generalized inversion must be used, i.e. 

A M. - S T [ SS T ] 1 AT . (7-14) 

Note that Eqs. (7-13) and (7-14) again assume the form of a linear 
guidance relation A u - P AT . 

Once Au has been determined, the updated controls are used to 
generate a new estimated trajectory to see if the target errors have 
decreased. This overall process is repeated until the target errors 
are made less than some specified tolerances, or until a maximum 
number of allowable iterations has occurred. In SIMSEP, the prin- 
cipal measure of target error is given by a so-called quadratic error 
function defined by 

* • £ K2M ’ 

where /^T(J) is the J*** component of the target error vector bnd 
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T tol (J) i# th* component of a vector of target error toler- 
ances. Convergence occurs In the nonlinear guidance whenever Q Is 
made less than one. 

Guidance divergence is said to have occurred if the quadratic 
error function is greater than the backup convergence criterion 
(A0K) after iterating a maximum number of times (NMAX). Divergence 
also occurs if the quadratic error function Increases c i three 
successive predicted corrections. As far as the Monte Carlo mission 
currently being executed, divergence is considered to be catastrophic, 
and the mission is ended. As a backup, weak convergence occurs if 
strong convergence falls to be satisfied but Q is less than A0K. 

In this way, the nonlinear guidance algorithm can be tolerant of 
"near misses" without bringing the mission to a halt. 

Another feature Included in the nonlinear guidance logic is 
the ability- to weight certain low thrust controls more than others. 
This permits the user an aided flexibility whereby he can effect a 
normalization of disparities in the units associated with controls. 

In addition, the user has the option to arbitrarily weight some 
controls more heavily, based upon knowledge and experience gained 
during the trajectory targeting process. 

Algorithms based on the computational steps outlined above are 
implemented in NLGUID. Both delta-velocity maneuvers and low thrust 
corrections are handled by essentially the same logic with A u in 
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4 1 

1 

I 

f 


c 



Eqs. (7-13) and (7-14) being a velocity update for impulsive guid- 
ance . 

7.4 Simulated Orbit Determination 

SIMSEP, in the strictest sense of the word, is not a complete 
"simulation" in that an explicit orbit determination process is 
not included in the computational algorithms. The problem of 
estimating a state vector is done by sampling a knowledge covariance 
in much the same way as it samples a control covariance or an ephem- 
eris error covariance. Simply stated, an augmented knowledge 
covariance is sampled to obtain an error in the estimated state 
vector, £, relative to the actual trajectory state, X^. Therefore, 
the estimated state vector is given by 

* h + * 

Likewise, parameters which have been augmented to the state and 
estimated during the orbit determination process are also computed. 

Typically, the knowledge error covariances which are read as 
input to SIMSEP have been computed in an equivalent G0DSEP run, 
they are equivalent in the sense that the same trajectory and 
sequence of guidance events are evaluated. With this procedqre, 
there is an added advantage of permitting a direct comparison of 
guidance results computed in G0DSEP and SIMSEP and their dependen- 
cies on the same state estimation results. 

There are several reasons why an explicit orbit determination 
capability has not been included in SIMSEP. Primarily, the 
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estimation process is not as subject to trajectory nonlinearities 

I 

as is the guidance process. This is because the estimation errors 
are generally small and well within the envelope of linearity. In 
addition, this method of simulating orbit determination minimizes 
the computational complexity of the program, while at the same time 
representing a cost effective means of performing an effective state 
estimation. 

7.5 Thrust Process Noise 

Digressing briefly to discuss the actual trajectory again, 
there is one very important process related to the generation of a 
real trajectory that is either ignored or modeled by an effective 
process in the other modes. This is the time -correlated thrust 
noise. These independent stochastic processes corrupt the commanded 
thrust controls, i.e., thrust magnitude, cone and clock angles, as 
small time-correlated perturbations, Each stochastic parameter is 
modeled in S1MSEP as a Gauss-Markov sequence which is computed during 
the actual trajectory integration. At time point tj^, the vector 
U, kn of stochastic parameters is given by 

-^-k+1 “ A + W 

— k+1 

where K, ^ has been evaluated at t^. t*. ^ is assumed to remain constant 

over the interval At = fc k+i -V with its effect being determined by 


the coefficient matrix, A 



110 


r - At/ t i 


0 


-At/T, 


0 


e * At ^‘* r N 




where 0*,, X>» ...» T are correlation times associated with 
i l n 

each corresponding stochastic parameter. is a vector of white 

noise terms which have statistics dictated by the requirement that 
the process remain stationary; namely 


O' 


Wj 


-2 At/Tj jr- 2 
(1 - e ) 


.th 


( r *■ 

where is the variance associated with the j v " component of 

the JL k+1 vector. During the integration, is evaluated at 

the start, the half-interval , and the end of a normal integration 
step. 

7.6 Guidance Execution Errors 

Once that a guidance correction has been formulated, the 
execution of that correction must be performed to affect the actual 

l 

trajectory. The commanded correction computed by the guidance is 
an Idealized set of control changes which are invariably cc rupted 
by execution errors. For a low thrust control change, the executions 
errors are actually built into the thrust process through the thrust 


i 



Ill 


biases and the dynamic process noise. However, for an. impulsive 
maneuver, an explicit set of logic to corrupt the ..coomanded deltc- 
velocity change must be implemented. 

In general there are three basic execution errors which are 
modeled for impulsive maneuvers: 1) pointing, 2) resolution, 

and 3) proportionality. Given the commanded delta-velocity vector, 

AV^ . in its heliocentric representation the in-and-out of the 
ecliptic plane angles are determined as, 

• C = tan’ 1 ( AV c (2) / AV c (3)) 

- sin’ 1 ( AV C (3) / IAV.I ). 

Specified pointing angle errors are sampled to formulate changes 
in the commanded angles (£'•< and ) to simulate orientation errors for 
the actual delta-velocity, -Mr Likewise, errors specified as 
proportional to the maneuver magnitude, £p , and as a minimum meas- 
urable resolution of a maneuver magnitude, 8 f , are also sampled and 
added to the comnanded correction. The actual velocity change is given as 

I A 7, y » | Av I + Sp + S r 


for the actual velocity magnitude. The actual velocity vector is 
given as 

6V a (1) - Uvj cos (<£ + £o< ) cos ($+S^> 

AV a < 2) " » & V sin(^+Sot) cos + 

AV 3) * 1 A V A 1 sin(£+Sp ) 
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APPENDIX 1 . 

% 

I 

9.1 Conic Equations For Position And Velocity In Elliptical 
And Hyperbolic Orbits 

Given: r , v , t and y 

— o -o’ o 

Find: £ and v at time t. 

From the initial conditions we can find the inverse semi-major 
axis 


1-1 

a r 


The mean angular motion 


L 

n m J— 


T®F 

and relationships for the eccentricity and eccentric anomaly for 
elliptical orbits (a>o) 

r • v 


e • sin E *===- 
o fTsT 


r v 4 

e • cos E “ — — — - i 

o p 


or for hyperbolic orbits (a<o) 

r • v 


e • sinh H - /- y a 
o 


r v * 

a • cosh H “ — — — 


-1 
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q;> 


To find the position and velocity along the conic orbit. Keplers 
equation must be solved for the change in eccentric anomaly. 

r 

Newton's method is used to solve the equation iteratively. Let 
Newton's method be given in the form 


f (* k ) 

*kH " x k ' F6T J . 


Then for elliptical orbits 


x - E - E 


( 


i’ 


f(x) ■ x + e*sin E (1 - cos x) - e*cos E sin x - n (t-t ) 
o o o 


and for hyperbolic orbits 

x * exp (H - H ) 


f (x) * H e-exp(H )*x + *se‘exp(-R )x+l - ln(x+i) - n(t-t ). 

o o o 

The position and velocity at time t in elliptical orbits is given by 
r -{l - f [l - cos(E-E o >] j r, + sin(E-E o ) 


- e*(sin E - sin E )r v 

o J -x 


- -liLlilsin (E-E ) • r + 
r r o -o 


{l - t [l - CO, <E-E o >1 J 


and for a hyperbolic orbit 


t 1 - f e l« sl > (H-V 1 ] } V n 

f e* (sinh) H-sinh H ) - *inh (H-H)l -v 

I O O J “V 

v . 8 inh (H-H ) r +fl - fTcosh (H-H) -l] [ 

— rr o -o L t w o *J 
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APPENDIX 2 

9.2 A Generalized 4th Order Runge-Kutta Algorithm With 
Runge'a Coefficients For A Matrix System Of First Order Differential 
Equations . 

The 4th Order Runge-Kutta formula for numerically integrating 
first order Differential Equations of the form 

y’ - f (x,y) 
is 

y k+l “ y k. + V (f 1 + 2 f 2 * 2 f 3 + V 

where h is the stepsize, x is the independent variable, y is the 
dependent variable and 

f i • f,( vV 

f 2 ' f + r • y k + -V- \ 

•V V f 2 > 

‘a' 1 ' (x k + r •» k + -V > 

'* ■ f ' (x i + V \ * W 

To generalize these equations for an mxn Matrix system of 
first order differential equations, write equation (1) as 

7 W (1 ,J) - y k (i,J) + ~ [f x (i» j) + 2-f 2 (i,j) + 2-f 3 (i,j) + 
f/i.j)] 
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and this can be written as 


where 


^ ■ Y k + r < F i + 2> f 2 4 *• f 3 + V - 


F i - F '<v V 


r 2 - F, (% + r- Y k + r- F i> 


? 3 - F, <*k + r • T k + r- f 2 > 


F 4 • T,i \+ \ F 3> 


> 
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APPENDIX 3 

I 

9.3 Newton's 3rd Order Divided Difference Interpolation 
Polynomial . 

Given: (x^ y^, (x 2> y 2 >, (x 3> y 3 ) and (x^, y 4 > 

Find: A third Order Polynomial that fits the given points 

Construct the following table 



where 


t X i+1' X J 


x i+r x i 


r x x x i . Ei+rvJ 

l X j+2 X j+1» X jJ x j+2 - Xj 


. i - 1,2,3 


. j • 1*2 


I *k+3 ,x k+2 .’St+l ,x kl 


1*11+3 * x k+2* x k+ll " t x k+2 ,x k l ,X k 1 




, k - 1 
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Newton's 3rd Order Divided Difference Interpolation Polynomial 
has the following form, 

* V 

5 (x) - yj+tx-x^- [x ,x ] + (x-j^) • (x-x 2 ) • [ x 3» x 2 » x il + 
(x-Xj) • (x-x 2 ) • (x-x 3 ) • [ x^ ,x 3 ,x 2 .x^ 
while a 3rd Order Polynomial can be written 


3 2 

y(x) ■ ax + bx + cx + d 

To find a,b,c and d, the coefficients of the x terms in the Newton's 

«* 

Divided Difference Polynomial can be equated to a,b,c and d, such 
that 



b * " (X3 1 X 2 , X^) i ^ Xj, x 2 , x L ] 

c - <Xy x 2 + x 3 , Xj.XjW - (Xj+Xj^)* [xj.Xj.x^ + [xj.x^ 
d - -(x lt x 2 ,x 3 )-a + x 3 x 2 [-j.Xj.x^ - [ x 2 ,x l^ + y 


Now a third order polynominal can be fitted to the four points and 
j can be determined from a given x or a maximum or minimum can 
be found from the following values of x. 


and 


-b +/b 2 - 3ac 

3a ^ 

_ -b -/b 2 - 3aji' 

3a 


is minimum 


is maximum 
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APPENDIX 4 

9.4 Analytic Expression! for Terms In the Fa Matrix 

In Section 4.5 it was shown that the augmented state transition matrix, 
♦a* is computed by integrating the matrix differential equation, 

*A * F A *A 

In order to efficiently integrate this expression, it is necessary to have 
analytic representations for the individual elements of F^. 

has been identified in equation 4*6 as a matrix of first order 
partial derivatives obtained by expanding the equations of motion. In 
concise symbolism, may be written as, 

F A - 

dM. 

where £a was defined in equation 4-5 and XA is the augmented dynamic state. 
For an analysis wh„re covariances are propagated by the state transition 
matrix, i.e., the STM mode, the (maximum) component vectors of X& are 
defined as: 


r 


s/c position vector 

V 


s/c velocity vector 



thrust control vector 

r 

HP 

m 

ephemeris body position vector 

V 

HP 


ephemerls body velocity vector 

Mv 


ephemeris body gravitational constant 

Ms 

» m 


solar gravitational constant 
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and tha corraapondlng matrix, in partitionad format, ia 


IW 

• 

0 

X 33 


0 


* 

0 

0" 


0 ” 


0 



» f 33 

0 

66 

*33. 

63 

„ k 33 

0 

mJ 

66 

- d 31- 

61 

-®31- 

61 

1 

0 

► 

0 1 

J 

36 | 

P 

1 ° . 

|33 

to 

0 

•] 

36 


[ 

I 31 

H 

bi 


0 

0 


o 


0 



0 


"0 



0 

0 

66 

0 

to m 

63 

_ p 33 

0 

66 


61 

_ S 31_ 

61 


0 

0 


r n 

0 


0 

0 ' 


0 


0 



0 

0 

26 

0 

23 

0 

0 

26 

0 


0 

_ J 


where the subscripts refer to the dimensions of the appropriate partitions, 


When the covariances are propagated by integrating the covariance differen- 
tial equation (see Section 4.6) in the PD0T mode, the augmented dynamic state 
vector is defined as 


r 


s/c position vector 

V 

m 

s/c velocity vector 

u 


constant thrust controls 

w 


time-varying thrust parameters (6x1) 

to m 


„ 4m 


and in this case is given as 



[° ° [°]b3 [ h ] 66 


(15x15) 
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In either SIM or PDOT node, the fully «ugt>9nted etate, as uaed by 
GODSEP, nay alao Include meaaurement parameters However, the terna 
appearing In corresponding to meaoureme'it parameters are zero because 
they do not affect the dynamic process. 

Specific matrices In are defined as follows: 
f^ * partlals of the s/c acceleration vector w.r.t. position 
components * d,a/d£ » 

^33 * partlals of the s/c acceleration vector w.r.t. the thrust 

controls * da/du , 

* partlals of the s/c acceleration vector w.r.t. the position 
of the ephemeris body * da/dr^ , 

djj ■ partlals of the s/c acceleration vector w.r.t. the gravitational 
constant of che ephemeris body * da/dfi » 

* partlals of the s/c acceleration vector w.r.t. the solar 
gravitational constant * da/dtf , 

p^j * partlals of the ephemeris body acceleration vector w.r.t. the 
position of the ephemeris body * » 

* partlals of the ephemeris body acceleration vector w.r.t. the 
ephemeris body's gravitational constant * da^/d^ , 

Sjj ■ partlals of the ephemeris body acceleration vector w.r.t. 
the solar gravitational constant • d a^ /d/y^ , 

n 36 * partlals of the s/c acceleration vector w.r.t. the time-varying 

thrust parameters ■ da/d w , and 

h,, * partlals of the time derivative of the time-varying thrust parameters 

w.r.t. the time-varying parameters * d£/dw . 





4 .. 


a 


122 


As noted from Its definition, the elements of the f 33 matrix are 
evaluated by differentiating components of the s/c acceleration vector, 
a., with respect to s/c coordinates, i.e.: 


33 


aa 

SI 


where £ is the sum of two contributing terms: the gravitational 
acceleration, £ , and the thrust acceleration, a^ . The partials of £ 
with respect to r are the components of the so-called gravity gradient 
matrix and are well-known, e.g. see Battin (Reference 4). In terms of 
s/c position vectors relative to the gravitating bodies, r^, the gravity 
gradient matrix is 


«/i>r - X [3 4 - r l * 33 ] "i /r i 


with the summation being performed for all bodies. In this equation, /t 

is the gravitational constant of the I*"* 1 body, and I^ * s a ‘*3x3) identity 
matrix. 

The second matrix comprising I s obtained by differentiating the thrust 

acceleration vector, as seen in the inertial frame, with respect to the 
heliocentric position of the s/c. Since the rotation matrix relating the 
body axis thrust vector to the inertial frame is dependent on the helio- 
centric position, it is necessary to differentiate components of the 
orthogonal rotation matrix as well as the components of the body axis 
thrust vector. Recalling from Section 4.1 (page 20) that the thrust 
acceleration is given by 

% = A , 
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this becomes upon differentiation, 

A 


*Jt .M a ' T) 

»r dr 


r a ! 

aj B 


ar 




ar 


®t 2 + a|!zi1 

r- J 


(A4-1) 


where £g, ^g and ^ are unit vectors defined by 


"b - I 1 |£| • 

% " i * ' I 1 “ -si 

A A 

= h*h' 


and 


A 

u 


r is, oi ccurse, the heliocentric position vector of the s/c, and Z is 

”S 


the unit vector of ecliptic direction cosines pointing toward the reference 
jtar for the cone/clock system. In terms of the unit vectors defined above 
and the unit vector for the reference star, the first three terms of the 
right hand side of equation A4-1 are 


&-f [ I33 -^ ** T ] 

A fA AT “1 

1_ Mb J B " *33 z 

dr ” R L J 

a^B = j [<M - K [c&bI 
ax [ax J j_dr J 


R is the magnitude of the vector cross product of _r and Is. The matrices 
Z, J, ana K are skew- symmetric matrices corresponding to the unit vectors 

A A 

Zs , jg and kg, respectively, and are defined as 





The last term in A4-1 reflects how the body axis components of the 
thrust acceleration vary with variations in the solar distance. Since the 
only quantity in each acceleration component which depends on the solar 
distance is the (scalar) power function, P(r), (See Page 16), the partlals 
of aL with respect to r, become 


I 



where the product of a^ 


1 

P(r) 
and 3P/<9r 


% mil 

ai 

I 

is a 3x3 matrix. 


The 3x3 g matrix is 

g 




- Stfabv 
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where A Is the transformation matrix from spacecraft cartesian to inertial 
coordinates (Section 4.1) and transforms thrust controls to spacecraft 

coordinates. 

_2i' « 

du 


(cone) 

-a’ 

y 

(cone) 

a* 

X 


0 


a^ i sin (clock) cos (cone) 

a’ -a* sin (cone) 

z 

for the cone/clock system, with a’ * thrust acceleration in spacecraft 


coordinates, and 

gj 1 = 


a’ 

x 


-a* 

y 


-a* sin£ cosV 

-a sin£ sinV 


a cos 


for the in /out of plane system. 
The 3x3 k matrix is 
k = ^ 

a 


I . 2£| [ 3 r r T -r 2 l1 - ^ 

X p V [ll P J R 

^ 3 R R T - R 2 I J 


Jis. * 

5 


where R ■ spacecraft position WRT the ephemeris body. 
The 3x1 d vector is 


t .p . .i 

dtip -s 


The 3x1 m vector is 


m 


dfta 
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u 


The 3x3 p matrix is 


d±2 = _ Up 


3r r 1 
“P “P 


r 2 I 
P 


The 3x1 q vector is 

„ . & . .3 


The 3x1 s vector is 




The 3x6 n matrix is 


r ■ r 


g i g 


The 6x6 h matrix is 




where •••» Tf, are P roc ess noise correlation times. 
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APPENDIX 5 


TOPSEP Inlection Modellt 


9.5.1 Injection Controls 

One set of controls which may be used in the targeting 
and optimization submode of TOPSEP is the initial state 5^* which 
defines a hyperbolic escape orbit relative to the launch planet, 
where 


2^ is the state of the S/C in cartesian elements immediately follow- 
ing injection from a parking orbit about the launch planet. Instead 
of using 5^ directly as control parameters, one may simulate the 
heliocentric injection process more realistically with a different 
set of injection parameters (See Figure 9-1, Page 128): 

1) r Q , the magnitude of the radius vector at the 
point of injection. 


2) i, 


the inclination of the parking orbit. 


3) Av, the magnitude of the velocity change from 

the parking orbit to a hyperbolic escape 
orbit. 

4) X , the angle locating the projection of Av 

vector in the parking orbit plane relative 


to r ♦ 



1 


*•**-■' , . . { I j I 

1 
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l 

5) 

♦ , the angle locating the Av vector out of 

v , l 


the parking orbit plane 


6) 

t Q , the heliocentric injection time. 


These Injection controls are augmented to the control profile in 
TOPSEP by correctly setting the following elements of the H array 
in the $T0PSEP namelist: H (9, 21), H (10, 21), H (1, 22), H (2, 22), 

H (3, 22), and H (4, 22). Refer to Page 15 in the User's Manual for 
additional information for implementation. The initial values of 
the injection parameters are not input in either the $TRAJ namelUt 
or the $T0PSEP namelist but are determined analytically from the 
initial state based on certain assumptions about the parking orbit. 

The injection state (r^, v q ) of the reference trajectory for 
any given iteration in TOPSEP is assumed to define a circular coplanar 
parking orbit (the eccentricity is very small but non-zero to accommo- 
date future program modifications for elliptical parking orbits). An 
in-plane Av is applied at time t which injects the S/C from the 
parking orbit to the hyperbolic orbit. Thus, the nominal or refer- 
ence parking orbit and injection controls are uniquely defined. 


r 

o 

i 

Av 

X 

* 


llol 

the inclination of the hyperbolic orbit 
J - Vpkj » where v^ k is the S/C's velocity 
in the parking orbit prior to injection 



0 

0 


I 

4 


a 
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The velocity in the circular parking orbit prior to injection is 


■Spk 




(r 

x v ) 

x r 

m 

— o 

-o 

— o 


l(r 

X V ) 

x r 


J 

-o 

-o 

- 


where A is the mass of the Earth. 

★ 

During th^ construction of the perturbed trajectories and 

trial trajectories in each iteration of TOPSEP, the injection controls 

are changed by some pre-determined amount so that a new cartesian 

state (r^ , v *) may be defined. When the injection controls r , i, 

— o — o o 

or t are modified^ new parking orbit must be computed. Changes to 
these injection controls affect both position r^ and velocity ; 
however, modifications to Av, X , or ♦ affect only velocity v^. 

If r , i, and t Q are changed to r Q , i, and t Q a new state 
(_/, v^ k ) may be found. First, the unitary angular momentum vector 
h, the node vector , and the longitude of the ascending node 12 

must be computed. If 


then 


and 


A 

h 




*hy 

hx 

0 


r xv, 
-o -pk 


h / h 


n ^ , 0* tan 1 (h u /-h y ) 


For the case in which the parking orbit is in the ecliptic, the node is 
defaulted to be the x axis so 


A 

n 




n« 0 radians . 


Finite differencing techniques are used to calculate injection control 
corrections for targeting purposes. 
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The angular position 0 of the S/C in the orbital plane as measured 


from the node is 


& traveled by changing 
injection time 


r / . 3 

< r o> 


( | 4- of S/C from node in j 

| nominal parking orbit j 

-1 / r . n \ 
cos | — o 1 


r - r 
— o o 


Finally, cosfl coS $ ~ sinfl sin 0 cos i # 

r # « r* sinfl cos A 4- cosQ sin 0 cos i 7 

-o o 

sin0 sin i' 

and 

V — — cosfl sin 0 + sinfl cos $ cos i 

fh sinfl sinfl - cosfl cosfl cos i* 

r ' a . ./ 

o -cosfl sin l 

is constructed by adding the delta-velocity vector A v* 

to the S/C*s parking orbit velocity v^* The vector Av^is computed 

in terms of in-orbit plane and out-of-orbit plane components. Let 
/ / / 

Av, X , and v|/ be modified injection controls. Then 
A v “ A v^ + A V£ + Av^ 

Av| “ A v 7 cos ( 'l' 7 ) cos ( X ) 

r b 

A Vj * Av 7 cos ( ty 7 ) sin ( X 7 ) h* x r^ 

|F^I 

' / t $ 

A v “ Av sin ( ♦ ) h 


where 


r ' x v ' . 
— o -pk 


^v 7 is the vector in the radial direction r^, Av^ * 8 c ^ e vector 
orthogonal to in the orbit plane, and Av^ is the vector normal to 
the parking orbit plane and in the direction of the angular momentum 


vector. Thus 


V + *■'- 


and a new initial state vector 


X* has been established. 


l 





Figure 9-1 illustrates how the injection parameters determine the 
new initial state vectors and vj. 



Figure 9-1. The new state vectors, r^ and v', as 

defined for perturbed and trial trajectories 
using injection controls. 


9.5.2 Tug Multiple-Impulse Orbit Transfer 


When the injection controls (r Q , i, t Q , Av, X , ) are 

applied in TOPSEP, the reference parking orbit is likely to change 
from iteration to iteration. In fact, it is possible that the park- 
ing orbit characteristic of the last iteration may not be attainable 
directly from an Earth based launch within realistic launch constraints. 
Therefore, it becomes necessary to consider the Interface between the 
SEP S/C and the launch vehicle in order to predict the "cost" of 
achieving the reference parking orbit. The cost may be estimated 
indirectly by determining the launch vehicle fuel budget to transfer 
the SEP S/C from some nominal inner parking orbit to the outer refer- 
ence parking orbit. If the inclination of the inner parking orbit is 
realistically constrained, an orbit plane change may be necessary to 
complete the transfer. As the angle of the plane change increases 
me estimated cost of the orbit transfer will increase dramatically. 

The estimated cost may then be used to distinguish between acceptable 
and unacceptable outer parking orbitr while simultaneously sizing the 
fuel expenditure for presently conceived or operational launch vehicles 
(or intermediate stages). An additional efficiency check may be made 
by comparing the fuel expenditure of the multiple Impulse orbit transfer 
with that of a single impulse injection from the inner parking orbit. 

The only requirement on the single impulse is that it provide the 

correct v vector. The single impulse calculation is not intended 
**• oo 

to be used in conjunction with the thrust control profile for targeting 
but rather as a standard for determining the inefficiencies of the 
multiple impulse injection process. 
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The launch vehicle simplistically modeled in TOPSEP is the 
* 

expendable space tug , Initially, the tug is in a circular inner 
parking orbit whose equatorial inclination is constrained due to 
bounds on the booster launch azimuth. The tug then performs the 
transfer to the outer parking orbit, the characteristics of which are 
described in Section 9.5.1. The outer parking orbit is also circular; 
however, it is assumed to have an inclination and ascending node equal 
to those of the hyperbolic escape orbit. The tug's path from the 
inner parking orbit to the outer parking orbit will be a coplanar 
Hohmann transfer if the required equatorial Inclination does not 
violate the launch azimuth constraints. Otherwise, the tug follows 
a modified Hohmann transfer (i.e., a regular Hohmann transfer in the 
inner orbit plane followed by a plane change and circularization at 
the line of intersection with the outer orbit plane). The equatorial 
inclination of this transfer orbit is either the maximum or minimum 
inclination bound such that the required plane change is a minimum. 

For the coplanar transfer the ascending node of the inner orbit is 
fixed; thus, the launch azimuth may be computed explicitly and tested 
for a constraint violation. For the plane change transfer the launch 
azimuth is fixed, and the inner parking orbit is uniquely determined 
when the minimum plane change condition is enforced. 

Figure 9-2 illustrates the selection of the inner orbit normal 
vector which minimizes the plane change for the transfer to the outer 
parking orbit.. 

*Although the launch vehicle is referred to as the space tug, in the 
remaining discussion, the specific vehicle is inconsequential to the 
sizing of the fuel budget (i.e., the vehicle may be a trans-stage. 
Burner II, Centaur, etc.) 
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ec 


normal vector of 
inner orbit mini- 



Figure 9-2 Selection of the Normal Vector for the Inner 
Orbit Which Minimizes the Plane Change 


The larger cone in the figure represents the locus of normal vectors for 
all possible inner orbits having the maximum equatorial inclination 
allowed by the launch azimuth constraints (AZMIN and AZMAX in $T0?SEP). 
The smaller cone is the locus of normal vectors for orbits having the 
minimum equatorial inclination, which in general will be equal to the 
latitude of the Kennedy Space Center (28.608 deg). If the normal vector 
of the outer parking orbit falls between the two cones, the parking 
orbits ere assumed coplanar. If the normal vector falls outside of this 
region, the inner parking orbit is characterized as follows: 
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(1) the Inner orbit is Inclined to the outer orbit by the 
angle 0 where 0 is the minimum angle between 

the normal vector of the outer crbit and the neat»»t 
cone, 

(2) the normal vector of the inner orbit becomes the 
projection of the outer orbit normal vector on the 
nearest cone. 

For both the Hohmenn and modified Hohmann transfers three 
distinguishable velocity Increments or Av's will occur. Figure 9-3 
illustrates the relative positions where these maneuvers are executed 
for the general case (r^ specified by RP1 in $TOPSEP). 



Figure 9-3 The Tug 3-Impulse Orbit Transfer 



The first impulse Av occurs at periapsis of the Hohmann transfer 

_ A 

at the line of intersection of the parking orbit planes. The second 
impulse Ay^ occurs at apoapsis of the transfer orbit and provides 
a velocity increment to circularize the orbit and ;o change orbit 
planes if necessary. At a point on the outer parking orbit prescribed 
by the injection position vector r^, a third impulse AY 0 executed 
which places the tug on a hyperbolic escape trajectory from Earth 
(Note that the injection impulse is exactly the same as the 

Ay vector discussed in Section 9.5.1). If both parking orbits are 
coplanar the line of intersection becomes ambigious so the impulse 
position vectors r ard r^ are assumed oriented such that 


r 

—a 


-r r 
a -t 


anJ. 




r 

o 


The total impulse at r^ becomes 

Av t = A^ + Av,; 

however, the impulses will always be referred to individually as in the 
orbital plane change example. 

The magnitudes of the velocity increments Av fl and Av^ for 
the general case are computed as follows. 



I 

1 


1 

I 

! 

I 
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S 



Since the first impulse does not initiate a plane change and since the 
velocity vectors before and after the impulse are aligned 



The second impulse may perform a plane change through the angle 0 . 
Thus, by applying the law of cosines 


Av, 


K 


) 2 + (v b -) 2 


- 2(v b ) (v fe ) cos 


♦) 


If 0 ■ 0 the above equation reduces to 
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The magnitude for the injection velocity increment Av q is 
dependent upon the state vector X^. The computations are discussed 
in detail in Section 9.5.1. 

Once Av^, an< * Av q ^ ave ^ een calculated the fuel 

budget for the space tug may be computed. The tug's dry weight (W ), 

tug 

the propellent weight an< * the tug's specific impulse (I gp ) 

are specified in the input namelist (TUGWT, TGFUEL, and TUGISP in 
$TOPSEP). The fuel required for the first maneuver (^£ ue ^) a is then 


^fuel^a 


(— (t?)) 


where 


= W + (W ,) + W 

tot tug fuel max sep 


W ge p is the weight of the SEP vehicle and g is the gravitational 
constant. It follows that the fuel for the second and third maneuvers 
are 



The required fuel budget to perform the orbit transfers including injection 
is then 


(W fuel ) 


tot 


- (W, ,) + 

fuel a 


(W r i>v + 

fuel b 


^fuel^o 
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If 




> 


tot 


®w 


max 


the outer reference parking orbit Is undesirable and measures must 
be taken to make the injection state and outer parking orbit realistic. 

Calculations for the single impulse injectior. from the inner 
orbit proceed as follows. The velocity at periapsis (v fl ) + necessary 
to obtain the hyperbolic excess veloc-ty is 

<V + * VmT ' ( Si> ' 

where a and e are the semi-major axis and eccentricity of the hyperbolic 
orbit. The parameters a and e are determined so that the vector v^ re- 
sulting from this orbit is the same as that obtained from the multiple- 
impulse injection process. 

The single velocity increment is then 



The required fuel expenditure is then 



A comparison of j) and (W will indicate the relative 

u a u tot 

efficiency of the multiple impulse injection process. 
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APPENDIX 6 

9.6 Control Weighting Schemes for TOPSEP 

Various weighting schemes are provided to allow the MAPSEP 
user flexibility in scaling controls in the TOPSEP mode (Chapter 5). 
The purpose of these schemes is to alter the contours of constant 
cost and constant target error in the control space so that the 
projected gradient algorithm may converge more readily. Convergence 
problems occur most often when elements of the control vector differ 
in units. For example, there may be considerable difficulty in 
finding a converged solution when thrusting angles (cone or clock) 
are selected as controls in addition to thrust phase times. Whereas 
the internal units for the angles and times are radians and seconds 
respectively, the corresponding elements of the sensitivity matrix 
may vary by several orders of magnitude. That is, the sensitivity 
of the targets to a change of one radian in thrusting angle is many 
orders of magnitude greater than the sensitivity of the targets to 
a change of one second in thrust phase duration. In the example 
just described one would find that the PGM algorithm would compute 
a control correction which would try to eliminate the target errors 
by large changes in the thrusting angles and very small changes in 
thrust duration. Unfortunately, large control changes are often 
invalid in the very nonlinear control space associated with low 
thrust trajectories. To alleviate this problem, the normalized 
control weighting scheme has been devised. The diagonal elements 
of the control weighting matrix are defined as 
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Application of this weighting matrix to the controls determines a 
sensitivity matrix whose elements are roughly the same order of 
magnitude. Thus, the removal of the target error is spread evenly 
among the selected controls rather than among Dnly a few* Hopefully 
all the cortrol changes will be small enough to be valid in the 
nonlinear control space. 

Another t'pe of convergence problem may occur. Sometimes 
elements of the sensitivity matrix vary by orders of magnitude 
even though the controls are all of the same units. For example, 
target parameters are much more sensitive to changes in thrusting 
angles early in the trajectory than they are to changes in these 
angles late in a trajectory. If the controls are not scaled, the 
PGM algorithm computes a control correction where changes in 
thrusting angles during early phases are unacceptably large and 
changes to Luc angles during later phases are undetectable. The 
following weighting schemes have been devised to spread the removal 
of target error more evenly among the selected controls, 
a) Sensitivity weighting 



b) Combined sensitivity, target error, and control weighting 


N 




Gradient weighting 



Averaged gradient and control weighting 




* 


+ 


0.1 

Uj + Gj ) 

1 2 
Gj 


where G is defined in c. 
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APPENDIX 7 


9.7 Integrated State Transition Matrices for Computing the 
Targeting Sensitivity Matrix 

Within the three basic modes of MAPSEP, trajectory guidance 
and/or retargeting represent one of the primary computational 
problems worked in the program. Whereas the logic controlling 
these calculations is, in general, straightforward and easily 
understood, the actual execution remains as one of the more costly 
computational operations to be performed. This is especially true 
in TOPSEP and SIMSEP where targeting over long trajectory arcs is 
done repeatly. In order to minimize computational expenses, an 
algorithm which uses state transition matrices integrated with 
the trajectory has been implemented and is used exclusively in 
GODSEP and SIMSEP for computing the targeting sensitivity matrix. 
In TOPSEP, the user has the option to use either this or the more 
expensive, but equivalent, numerical differencing algorithm. 

The targeting sensitivity matrix, S $ is a matrix of linear 
partials relating variations in the control variables, £ u, to 
variations in the targets, A T , according to 

A 1 * s A e • 

I 

Looking at the targeting sensitivity matrix in more detail, it is 
seen to be of the form, 

£ ( T l» T T n ) 

S ! - . -■■■ I 

v' ( V u 2 V 


where the T f s are selected target variables and the u f s are controls. 



Typical target variables Include X^, Y^, Z { , B*T, B*R, etc. which are 
all evaluated about the final trajectory state. Typical low thrust 
control variables are the thrust phase stop time, thruster throttling, 
cone angle, angle in the same or different thrust control 
phases. 

To provide a conceptual understanding of how S is evaluated 
from trajectory information generated in the augmented state transi- 
tion matrix, it is convenient to consider the trajectory segment 
depicted below where time points k, 



k + 1, .... and f are shown bounding thrust control phases n, n + 1, 
etc. In the figure, time point f denotes the target time and represents 
the trajectory stopping condition. Considering a specific thrust con- 
trol phase, say n, it is recalled from Section 4-5 that the augmented 
state transition matrix generates partials of state variations of time 
k + 1 with respect to state changes at k and control variable changes 
interior to thrust phase n. In particular, these partials are contained 
in the | and § u partitions cf the augmented matrix and may be sym- 


bolically written as, 



I 

I 
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A<«. y. *■ v x . v v . y,) k + i , 

<S (*, y, X. v x , X y , v 2 ) k 


and 


0 

u 


^(x, y, z, v x> v v> v z ) k ± 1 . 
«^(u r u 2 , ~ t u 4 ) 


The u's correspond to the phase stop time, throttling, cone angle, 

t h 

and clock angle for the n phase. (Cone and clock angle rates are 
excluded from this treatment since their partials must be obtained 
by numerical differencing.) 

If u^ (for example) in thrust control phase n is specified as 
an active control for the targeting event being considered, then the 
action of on the state vector at f is computed by pre-multiplyinj^ 
the appropriate column from 0^(n) f° r phase n by all intervening 
| 1 s, i.e. 


Ik 

^u 2 (n) 


f f, k+3 


‘k+3, k+2 Ik+2, k-H 


(n) 


Hence, an augmented 0 u> say 8^, can be constructed by storing and 
pre-multiplying appropriate columns from the 0^'s 88 they are computed 
during trajectory integration. The resulting gives the partials of 
the final state with respect to the active controls occurring in the 
various thrust phases between k and f and may be written as 


A 

§ 

u 


^ (x, y, z, v x , v y , v z ) f 

<* ( V u 2 V 


i 



t 

1 


i 


I 
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The final component necessary to compute the requisite S matrix 
Is the evaluation of target variable partlals with respect to the 
final state, This Is most expediently done by a numerical dif- 

ferencing algorithm to generate the differential point transformation 
matrix, . The matrix can be written as 


/v\ , — ^ Tl> T JlL Isi 

\ z » V x» v V z >f 


Now S is seen to be 


V- 


In TOPSEF where initial state conditions are also permitted as 
control variables, it is necessary to extend the above procedure but 
not the computational method. For this special case, only the chained 
$ ' s are needed to compute the partlals of final state with respect 
to charges in the state at k. Clearly, selected columns from §, , 

•tjK 

which is defined by 

If,k * }f,k+3 Ik+3, k+2 Ik+2, k+1 Ik+1, k 


may be augmented to the previously defined to give the requisite 
targeting sensitivity matrix to be used by TOPSEP. 
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